In this paper, we report on a series of mobile‐bed laboratory flume experiments designed to quantify the effect of sediment transport rate on the development of bed surface structure in gravel‐bed rivers. Most of the studies referred to above concentrate on describing the structural properties of equilibrium channel beds and as a result, they sample the bed relatively infrequently. A novel aspect of this study is the adoption of a high sampling frequency in an explicit attempt to capture the evolution of gravel‐bed structure as equilibrium bed conditions develop. Another innovation is the use of a regional‐residual separation technique to isolate different components of the bed surface topography; specifically the contributions made at mesoscale and grain scale. The primary goals of the paper are to: (1) describe and account for the structural characteristics of mobile armors formed in a laboratory flume; (2) quantify the evolution of mobile armor structure from an initial, screeded, bed surface condition; and (3) explore how mobile armor structure is influenced by transport rate.

Differences between the structural characteristics of screeded and water‐worked gravel beds have been attributed to the development of streambed armors [ Nikora et al ., 1998 ; Cooper and Tait , 2009 ]. Several workers have reported the results of flume experiments showing the simultaneous development of bed structures and static armors (those formed under conditions of zero sediment flux; Church et al . [ 1998 ]). Under these conditions, the structural properties of static armors have been related to the armoring discharge [ Aberle and Nikora , 2006 ] and the proportion of sand in the bulk sediment mixture [ Curran and Waters , 2014 ]. Other workers have considered the structural properties of mobile armors formed under a variety of upstream sediment supply conditions [e.g., Hassan and Church , 2000 ; Marion et al ., 2003 ] and have highlighted some differences with the static armor case. Mao et al . [ 2011 ], for example, found that static and mobile armors differed in their response to increasing discharge: static armors displayed limited variation in vertical roughness scale and became less topographically complex whereas mobile armors exhibited increased roughness length scales, less organization, and larger cluster structures. The differences were attributed to the active transport and greater protrusion of coarse bed material in the mobile armor case. The literature also suggests that bed surface structure should be influenced by sediment sorting. Reid et al . [ 1992 ] and Wittenberg [ 2002 ], for example, report field results that demonstrate that the occurrence of pebble clusters increases as bed material sorting decreases (i.e., as sorting indices increase), and Church et al . [ 1998 ] note that bed structures are rarely reported in flume studies using narrowly graded sediments. This is counter to Strom et al . [ 2004 ], who recorded the development of pebble clusters in uniform sediments in their flume experiments. Other workers have stressed hydrological controls. Wittenberg [ 2002 ] found that the density of pebble clusters was lower in dryland ephemeral rivers than in humid‐temperate perennial rivers and attributed the difference to the flashy discharge regime and absence of base flow in dryland channels, which restricts the period for cluster development. The importance of flow regime is supported by the work of Haynes and Pender [ 2007 ] and Ockelford and Haynes [ 2013 ], who demonstrate significant changes to packing arrangements and particle orientations during the periods of subthreshold flows between competent flow events, and by Mao [ 2012 ], who documented time‐dependent variations in bed structure during the passage of stepped hydrographs. Other studies suggest that a temporal dimension to bed structural development may exist. Increases in bed roughness with flow duration have been attributed to structure development, as well as bed surface grain‐size changes, by Pender et al . [ 2001 ], Marion et al . [ 2003 ], Mao et al . [ 2011 ], and Ockelford and Haynes [ 2013 ], and Marion et al . [ 2003 ] have further shown that different bed structures can form over different timescales.

Conventionally, the surfaces of gravel‐bed rivers are described and characterized by particle size and, as a result, there is a large literature concerning sampling techniques and standards for determining particle size characteristics of streams and rivers [e.g., Bunte and Abt , 2001 ]. However, river bed roughness and stability cannot simply be attributed to particle size because the development of bed structure—variations in the spacing, packing, and geometrical arrangement of particles [ Laronne and Carson , 1976 ]—evolves concurrently with armoring as particles in traction come to rest in stable positions, often against other stable particles. Recognizable particle configurations tend to develop, reflecting the arrangement of individual particles, for example, alignment or imbrication [ Johansson , 1976 ; Qin et al ., 2012 ], and the organization of particles into coherent bed forms such as pebble clusters [ Brayshaw , 1984 ] and stone cells [ Church et al ., 1998 ]. There is a growing understanding of how different bed structures contribute to bed roughness and stability through the generation of turbulence [ Lacey and Roy , 2008 ; Papanicolaou et al ., 2001 ; Hardy et al ., 2010 ; Tan and Curran , 2012 ; Curran and Tan , 2014 ], the promotion of form drag [ Hassan and Reid , 1990 ; Clifford et al ., 1992 ; Hassan and Church , 2000 ], modifications to particle exposure and grain pivot angles [ Kirchner et al ., 1990 ; Hodge et al ., 2013 ], and the spatial distribution of cluster forms [ Piedra et al ., 2012 ]. However, relatively little is known about the factors that are conducive to the development of bed structures in gravel‐bed rivers and whether different conditions promote the development of different types or degrees of bed structuring.

The bed surface of an alluvial river is the primary interface between the flow of water and the sediment available for transport. As such, it is a critical component of the fluvial system, exerting a fundamental influence on near‐bed hydraulics [ Hardy et al ., 2010 ], resistance to flow [ Powell , 2014 ], and sediment transport rate and grain‐size distribution [ Gomez , 1995 ]. Moreover, as a result of feedbacks between the bed, the flow, and the transported sediment, the channel bed is free to adjust to changes in discharge [ Parker et al ., 2003 ] and sediment supply [ Dietrich et al ., 1989 ]. It thus represents a potential degree of freedom in river self‐organization and adjustment [ Ferguson , 2008 ; Ferguson et al ., 2015 ].

Directional variograms of gravel‐bed river roughness often have a characteristic form that can be divided into three regions [ Nikora et al ., 1998 ; Nikora and Walsh , 2004 ; Cooper and Tait , 2009 ; Mao et al ., 2011 ]. At large spatial lags there is a saturation region where ( h x , h y ) → s 2 , and the data exhibit no spatial dependence. At short spatial lags there is a scaling region where the variogram takes the form of a power law γ = ch 2 H , where c is a constant and H is the Hurst scaling exponent. In between these two regions is a curved, transition region. Taken together, the scaling and transition regions define the lag distance termed the range ( h r ) of the variogram. The spatial lags at the upper ( h 1 ) and lower ( h 2 ) limits of the scaling and saturation regions, respectively, are important properties of the variogram: the former defines a characteristic, horizontal, bed‐roughness scale while the latter defines the correlation length of the bed surface elevations [ Cooper and Tait , 2009 ]. The value of the Hurst exponent ( H ) is also significant because it provides a scale independent measure of positive autocorrelation (or persistence) of bed elevations for 1 ≤ h ≤ h 1 mm [ Robert , 1988 ]. Following Bergeron [ 1996 ], objective estimates of these parameters were defined by linear spline interpolation [ Smith , 1979 ] whereby h 1 and h 2 are defined by the two knots of a three degree spline and the scaling exponent H is defined by the gradient of the line fitted to the scaling region (1 ≤ h ≤ h 1 ). The model was implemented in log‐log space using the PROC NLIN function of the SAS software package [ Freund and Littell , 2000 ]. The slope of the third line segment (fitted to the saturation region; h > h 2 ) was constrained to be zero because the semivariance in this region approximates the sample variance.

Variogram analysis [e.g., Oliver and Webster , 1986 ; Robert , 1988 ; Butler et al ., 2001 ] was used to provide a further view on topographic variation with scale. In this study, the spatial dependence of the relief across the experimental surfaces was evaluated using directional variograms constructed in the downstream ( γ ( h x , 0)) and cross‐stream ( γ (0, h y )) directions and using 2‐D variogram surfaces ( γ ( h x , h y )), where γ is the semivariance calculated for a number of spatial lag distances ( h ) measured in the downstream ( h x ) and cross‐stream ( h y ) directions (equation (7) in Table 2 ) for 1 ≤ h x ≤ 490 mm and 1 ≤ h y ≤ 190 mm. The upper limits of h x and h y approximate half of the downstream and cross‐stream dimensions of the surface, respectively, and ensure that sufficient pairs of points can be sampled from the surface at the longest lags in order to calculate γ accurately [ Klinkenberg , 1994 ]. The lower limit of h was chosen to minimize the influence of the subgrain‐scale microtopography [cf. Butler et al ., 2001 ].

Probability density functions (PDFs) of surface elevations and their properties (i.e., range ( R ), standard deviation ( s ), skewness ( Sk ), and kurtosis ( Ku *)) were calculated using equations (1)–(4) in Table 2 . Although such parameters are relatively crude descriptors of bed microtopography, several workers have shown them to reflect surface‐forming processes [e.g., Marion et al ., 2003 ; Aberle and Nikora , 2006 ; Mao et al ., 2011 ; Mao , 2012 ]. The area ratio ( A R ) was calculated using equation (5) to provide a measure of surface roughness [ Hobson , 1972 ; Grohmann et al ., 2011 ; Brasington et al ., 2012 ]. This is a ratio between the actual and plan surface areas [cf. Grohmann , 2004 , Figure 10] and is thus large for rough surfaces and close to one for smooth ones. Equation (5) represents a standard method of calculation based on a surface derivative [e.g., Horn , 1981 ; Grohmann , 2004 ] and is equivalent to, but simpler to implement than, that of Brasington et al . [ 2012 ] and Grohmann [ 2004 ]. The inclination index ( I ) of Smart et al . [ 2004 ] was calculated to highlight the development of grain imbrication (equation (6) in Table 2 ). The index reflects the balance between the number of downstream and upstream facing slopes and is negative for an imbricated bed [see Hodge et al ., 2009b , Figure 2].

The topography of the test section was captured using a hand‐held laser scanner (Faro Laser Line Probe®) mounted on a seven‐axis articulated arm (Faro ScanArm®). Although the arm provides for the free and easy movement of the laser over and around the complex geometry of the bed, it proved difficult to use within the confines of the flume walls, hence the need for the removable test section. Once the test section had been removed from the flume, the surface was scanned in numerous overlapping short strips up to 60 mm wide with a point density of up 640 points per scan width and precision of 0.086 mm. The resultant point clouds were aligned and registered to a local coordinate system as a single point cloud using Polyworks® proprietary software. Although the resultant point clouds are in 3‐D, algorithms that produce fully 3‐D surfaces from 3‐D point clouds are complex and have significant processing requirements. Following Hodge et al . [ 2009a ], therefore, the point clouds were filtered by removing points lying directly beneath, or close to, a high point with the same ( x , y ) coordinates using a point spacing of 0.05 mm. In effect, the filter ensures that each ( x , y ) location in the point cloud becomes associated with a single elevation ( z ). Digital Elevation Models (DEMs; hereafter referred to as the source DEMs) were then created by interpolating the filtered data onto a 0.15 mm grid using the inverse distance weighting algorithm in ArcGIS®. This resolution provides multiple points on any of the smallest grains on the surface, thereby maximizing the chances of including them in the bed topography. The source DEMs were cropped to a rectangular footprint of 990 × 390 mm and then linearly detrended to remove any large‐scale slope. These DEMs are subsequently referred to as the measured DEMs.

Bed surface grain‐size distributions were determined from photographs of the removable test section each time it was removed from the flume for scanning. The surface grain‐size distribution was sampled by superimposing a grid over the photographs and noting the colors of grains falling under the grid intersections [ Wilcock and McArdell , 1993 ]. Each sample comprised 400 point counts taken from the same 400 × 400 mm area of the test section using a 20 mm square grid.

For each experiment, the flume was set to the desired slope and the initial surface prepared by running the flume for 30 min at 10 mm flow depth (0.003 ≤ τ * ≤ 0.007) in order to remove any artifacts of the screeding process (e.g., any particularly unstable grains) and to displace any air. The bed was drained, and the test section was removed from the flume. The surface topography was then captured and an area in the center of the test section photographed for the determination of surface grain size. The test section was then replaced, and the flow restarted. This procedure was repeated a further eight times during each experiment producing photographs of the surface grain‐size distribution and surveys of the surface topography at 0, 30, 60, 120, 180, 360, 620, 900, and 1200 min. These intervals were chosen on the basis of pilot runs in which the transport rate initially declined very rapidly, suggesting that the greatest change in surface conditions would happen at the beginning of the experiment [e.g., Pender et al ., 2001 ]. Bed load data were collected in association with the sediment recirculation protocol discussed above. That is, prior to recirculation, sediment from the trap was weighed to give estimates of bed load transport rate at t = 10, 20, 30, 45, 60, 120, 180, 360, 620, 900, and 1200 min. For each experiment, the average value of bed shear stress was calculated as τ = ρgRS , where ρ is density of water, g is acceleration due to gravity, R is hydraulic radius, and S is slope (corrected for sidewall effects using the method of Williams [ 1970 ]). For each surface, values of dimensionless shear stress were calculated as τ * = τ /( ρ s − ρ ) gD s 50 , where ρ s is the density of the sediment (taken to be 2650 kg m −3 ) and D s 50 is the median particle size of the bed surface.

Sediment collected in a sediment trap at the downstream end of the working section was periodically reintroduced at the upstream end of the flume to provide a form of sediment recirculation appropriate for coarse‐grained materials [ Wilcock and McArdell , 1993 ; Mao et al ., 2011 ]. For the first 60 min of each experiment, transport rates were high due to the unconsolidated nature of the bed. Accumulated sediment was therefore reintroduced more frequently during these early stages of the experiments. Reintroductions took place every 10 min for the first half an hour, every 15 min during the following half an hour, and thereafter at t = 120, 180, 360, 620, and 900 min. Pilot runs demonstrated that during the first 120 min of each experiment, the intervals for bed load reintroduction were long relative to rates of bed adjustment. In this case, reintroduced sediment was added to a system which had already adjusted to a lower transport rate than that which moved it into the trap and, instead of maintaining an equilibrium, created a static deposit at the upstream feed area. Reintroducing the full volume of sediment produced in the previous interval was therefore inappropriate. Instead, we used pilot runs to establish the amount of sediment that was anticipated to be produced in each interval. Then at each reintroduction period during the first 120 min the amount of sediment that was reintroduced was equivalent to the yield anticipated for the next sampling period (care was taken to ensure that the grain‐size distribution of the reintroduced sediment was representative of the grain‐size distribution of the sediment caught in the trap).

Flow depth ( Y ) was set at 120 mm throughout each experiment to maintain a width:depth ratio of 5, thereby limiting the development of secondary flow cells [ McLelland et al ., 1999 ]. The relative roughness ( Y/D 84 ) of the experimental mixture was 5.2 and, therefore, at the lower end of values that are typically observed in gravel‐bed rivers ( Y/D b 84 > 5; [ Bathurst , 1993 ]). Experiments were conducted at three slopes ( S = 0.006, 0.009, and 0.015) which generated sidewall‐corrected, dimensionless shear stress values ( τ * ) for the initial ( t = 0 min) sediment surfaces of between 0.031 and 0.088 (Table 1 ; 1.0 ≤ τ * / τ * c ≤ 2.93; sidewall correction after Williams [ 1970 ]). Equivalent values for the final surfaces ( t = 1200 min; Table 1 ) were slightly lower as a result of surface coarsening and ranged from 0.023 to 0.047 (0.77 ≤ τ * / τ * c ≤ 1.57). Each experiment lasted 1200 min. Since pilot runs had established that fluctuations in transport rate decreased markedly over the first 360 min of each experimental run and had fallen to ~1% of the initial fluctuations after 600–900 min, this ensured that the final surfaces generated in each experiment were produced by a sustained period of quasi‐steady transport, [see also Mao et al ., 2011 ]. In all cases, flows were steady and uniform through the test section with a relatively stable water surface.

A rigid tray (1000 × 400 mm), formed from 8 mm thick stainless steel, was buried in the test sediments with its upstream edge 3 m from the flume inlet. The tray allowed removal of a portion of the bed, the test section, to facilitate measurements of surface topography, which were made outside the flume. Great care was taken when draining and refilling the flume and removing and replacing the test section so that disturbance to the bed surface was minimized, both inside and outside of the tray. At each sample interval, the pump was stopped and the flow was allowed to drain naturally. This was observed to stop transport immediately without causing any obvious disturbance to the bed surface, e.g., by drawing fines down through the bed surface. The test section was not removed from the flume until it had entirely drained so as to further minimize disturbance. The tray was removed using an overhead lifting device that enabled the tray to be lowered onto a bench beside the flume for scanning. Insertion of the tray back into the bed was aided by a retainer of mesh fencing that prevented collapse of the surrounding bed when the tray was removed. Measurement of bed surfaces in this manner has successfully been used before [ Marion et al ., 2003 ; Cooper et al ., 2008 ; Ockelford and Haynes , 2013 ]. Refilling the flume was undertaken carefully with the bed flooded slowly before discharge was increased to the required level.

The sediment mixture was made up from natural quartz‐density river sands and gravels ranging between 1 and 45 mm in diameter. The sediment was sieved into 1/2 phi size fractions and recombined into a grain‐size distribution with an inclusive graphic sorting coefficient ( σ b ) of 1.5 [ Folk and Ward , 1957 ] and D b 50 of 8 mm (Figure 1 ). The size distribution of the sediment mixture was not directly scaled from field samples, but the sorting value is representative of the lower limit of values found in 148 subsurface grain‐size distributions from nine gravel‐bed rivers, including Fraser, Sukunka, and Peace Rivers and Carnation Creek (Canada); Colorado River and Redwood Creek (USA); Drôme River (France), Ngaruroro River (New Zealand); and the River Wharfe, UK. These grain‐size distributions were collected by us or were provided by colleagues (see acknowledgments). For these samples from predominantly large, piedmont, rivers with a riffle‐pool style, the median sorting coefficient was 2.58 and the fifth percentile was 1.56, such that the degree of sorting in the sediment mixture used here falls in the lower tail of the natural distribution (since the proportion of the field grain‐size distributions < 1 and > 45 mm were, on average, similar (cf. 0.15), the experimental mixture was not truncated more at one end of the grain‐size distribution than the other). The well‐sorted sediment helps to isolate sediment transport rate as the dominant treatment factor in the experiments, given the potential increase in structural development as sorting decreases [e.g., Reid et al ., 1992 ]. Following Wilcock and McArdell [ 1993 ], all the grains in each size fraction were painted a specific color so that surface grain‐size distributions could be determined by counting the numbers of differently colored particles exposed on photographs of the bed.

Experiments were performed within a glass‐sided, 8.2 m long, flow‐recirculating flume of rectangular cross section (0.6 m × 0.5 m). To prevent scour and to induce turbulent boundary conditions at the flume entrance, 2.1 m of immobile sediment was placed directly downstream of the water inlet. No measurements were made in the final 1.7 m of the flume in order to avoid drawdown effects. The flume, therefore, had an effective working length of 4.4 m which was covered with thoroughly mixed mobile test sediment to a depth of 90 mm (equivalent to ~ 11 D b 50 and ~ 4 D b 84 , where D bx is the x th percentile of the bulk experimental sediment mixture) which were screeded flat so the sediment surface was parallel to the flume floor.

Three flume experiments (e1–e3) were conducted with a sediment recirculation protocol to compare the temporal evolution and surface properties of mobile armors formed at different flow strengths. Each experiment was repeated once (e1r–e3r) to assess the replicability of the experiments. The flows were designed using the concept of transport intensity, τ * / τ * c , where τ * is the dimensionless shear stress and τ * c the value of τ * at the threshold of motion, which is typically low in most gravel‐bed rivers, even at bank full flow ( τ * / τ * c ≤ 1.6; [ Parker , 2006 ]). In this study, different flow strengths were generated by keeping the flow depth constant and varying the flume slope. We selected the slope of each experiment to generate a range of representative transport intensities under the assumption that τ * c = 0.03, a value that is toward the lower end of estimates published in the literature [ Buffington and Montgomery , 1997 ; Petit et al ., 2015 ] but typical of bed surface conditions that develop in laboratory flumes at low τ * [ Ferguson , 2012 ]. The experimental sediment comprised a mixture of fluvial sands and gravels. The development of the bed surface condition was evaluated at the beginning and end of each experiment and at seven time points during the experiments using measurements of surface grain size and topography. Bed load transport rates were also measured by periodically sampling the sediment exiting the flume. Table 1 provides a summary of the test conditions for the three treatments that made up the experimental program. Note that no data on the final bed surface condition is available for experiment e3r and that this experimental run, therefore, is truncated at the preceding sample interval ( t = 900 min where t is the elapsed time of the experiment).

3 Results

3.1 The Bed Load Dynamic The character of the bed load dynamic in each of the experiments is shown in Figure 2. The different experiments share the same pattern of variation in that bed load transport rates were initially high, but decreased rapidly during the early stages of each run, with much of the adjustment occurring in the first 120 min (Figure 2a). Some further adjustment occurred over the next 240 min as transport rates declined less rapidly; thereafter, transport rates averaged over 260–300 min stabilized (the average difference in transport rate between sampling intervals is −0.12 g m−1 s−1), though some experiments (most notably experiment e2) continued to decline indicating that steady state conditions were approximated. The final transport rates varied by a factor of 7 between the three experiments and, as expected, increased with imposed shear stress. The increase in transport rate between experiments e2, e2r and e3, e3r (where τ* ≈ 0.028 and 0.046, respectively) was not, however, as great as that between experiments e1, e1r and e2, e2r (where τ* ≈ 0.023 and 0.028, respectively). The reduced sensitivity of transport rate at the higher flow strength is unexpected and may reflect either a reduction in transport stage (τ*/τ* c ) as a result of a structure‐induced increase in entrainment thresholds and/or a roughness‐induced reduction in the available shear stress or a limitation of sediment supply as a consequence of the manner by which sediment was recirculated. The latter, however, seems more likely given that our analyses of the bed surface condition described below indicate little difference in the surface structure and roughness of the beds developed in experiments e2 and e3. Changes in the bed surface grain‐size are shown in Figure 2b. During the experiments, the bed surface developed a coarse surface armor layer as the proportion of fine and coarse grains in the surface size distribution decreased and increased, respectively. As with the transport rates, most of the adjustment in surface grain size occurred in the first 120 min and had largely been completed by t = 360 min. Figure 2c compares the D s50 and D s84 of the final armored beds. It suggests that the armor layers coarsened from e1 to e2 but that little further coarsening was observed at the higher shear stress applied in e3. Inspection of the bed load grain‐size distributions suggests that there was little preferential winnowing of fines from the surface which, under the sediment recirculation protocol, might have promoted the coarsening of the surface. Instead, surface coarsening is attributed to a gravity‐driven kinematic sorting process in which finer grains preferentially move into voids left behind by entrained grains [Parker and Klingeman, 1982]. Since the kinematic sorting process is not thought to depend explicitly on the rate of transport, mobile armors that are invariant with flow strength are the expected pattern of equilibrium bed surface adjustment in a laboratory flume when the sediment output is recirculated [Wilcock and McArdell, 1993; Wilcock, 2001; Mao et al., 2011]. Figure 2 Open in figure viewer PowerPoint D 84 of the final bed surfaces with shear stress for each experiment and its replicate. In Figures D 50 and D 84 of the six initial experimental surfaces; the solid lines represent the D 50 and D 84 of the final surfaces when averaged for each experiment and its replicate. (a) Variation in bed load transport rate and (b) median particle size of the bed surface (key as in Figure 2 a) with time and (c) variation in median particle size andof the final bed surfaces with shear stress for each experiment and its replicate. In Figures 2 a and 2 b, the data points are shown for one of the experiments (e1) for reference. In Figure 2 c, the dashed horizontal lines represent the averageandof the six initial experimental surfaces; the solid lines represent theandof the final surfaces when averaged for each experiment and its replicate.

3.2 Visualization of the DEMs Some examples of the measured DEMs of the experimental surfaces are shown in Figure 3. The screeded beds (Figures 3a–3c) had a flat uniform surface with random grain‐scale variations in bed elevation about the mean. The impression was of a loose arrangement of particles. In contrast, the water‐worked surfaces comprised well‐packed, imbricated particles and displayed considerable heterogeneity at the mesoscale (10–100 mm) with patches of lower and higher bed elevations clearly discernible (Figures 3d–3i). The water‐worked topography developed rapidly with the onset of sediment transport (it was clearly evident in the surfaces at t = 30 min (Figures 3d–3f)) and the mesoscale topography was more pronounced in experiments e2 and e3 which were run at higher shear stresses (Figures 3g–3i). The spatial structure of the mesoscale topography varied over the duration of each experiment and was reminiscent of the longitudinal ridges and troughs which appear to characterize many water‐worked beds [Marion et al., 2003; Hodge et al., 2009b; Cooper and Tait, 2009; Bertin and Freidrich, 2014] and may reflect the development of secondary flow cells [McLelland, 2013]. At the grain scale, clusters of particles were observed in the water‐worked beds, but they did not appear to be a significant feature of the bed topography. This may reflect the narrow grading of the experimental mixture. Figure 3 Open in figure viewer PowerPoint Digital elevation models of bed surface topography for experiments e1, e2, and e3 at (a–c) t = 0, (d–f) 30, and (g–i) 1200 min. The bed area is 990 × 390 mm, and the flow direction is from left to right.

3.3 Probability Density Functions of Bed Elevations and Their Moments During each experiment, the range and standard deviation of bed elevations increased, the number of observations around zero mean decreased, and the PDFs become positively skewed with broader, flatter peaks. The dynamic of this adjustment is illustrated in Figure 4a which shows the temporal evolution of bed elevation PDFs in experiment e2. The initial, screeded bed was characterized by a narrow, peaked distribution with slight negative skew and positive kurtosis. These characteristics were shared by all of the start‐up beds and indicate consistency in the preparation of the initial surfaces. The most significant adjustment in the shape of the PDFs occurred in the first 30 min as the range and standard deviation increased by 7.4 mm and 2.2 mm, respectively, and the maximum probability density decreased from 0.11 at zero mean elevation to 0.06 at an elevation of −2 mm. Thereafter, changes in the shape of the PDFs are relatively modest. This pattern of adjustment was observed in all experiments, though the magnitude of the adjustment varied between treatments. A comparison of the PDFs for the final surface of each experiment indicates that higher flows generated broader, flatter distributions with heavier tails (Figure 4b). Figure 4 Open in figure viewer PowerPoint t = 0 and 1200 min for reference and to highlight the initial and final surfaces. Probability density functions of bed surface elevations for (a) all the surfaces of experiment e2 and (b) the final surface of each experiment. In Figure 4 a, the data points are shown for= 0 and 1200 min for reference and to highlight the initial and final surfaces. Further insight into the nature of the bed surface adjustment under different flow strengths can be obtained by considering the evolution of the statistical moments of the bed elevation distributions. These are generally consistent between replicates of each experiment. For each experiment, Figure 5 shows how standard deviation, skewness, and kurtosis varied over time (Figures 5a–5c) and how the moments calculated for the final, equilibrium surfaces varied with dimensionless shear stress (Figures 5d–5f). Much of the adjustment in standard deviation occurred in the first 120 min of each experiment with little change occurring thereafter (Figure 5a). The standard deviation of the equilibrium surfaces increased from 5.4 mm at the lowest flow to 7.1–7.7 mm at the intermediate flow, with only a modest further increase at the highest flows (s ≈ 8 mm; Figure 5d). Figure 5b shows that, with one exception (e2r), the bed elevation PDFs shifted from negative to positive skew within the first 30 min of each experiment, with surfaces exhibiting either near‐symmetrical (Sk ≈ 0) or positively skewed PDFs thereafter. Kurtosis also adjusts rapidly during the early stages in the experiments, reflecting a shift to more platykurtic distributions (Figure 5c), although the variation thereafter is more pronounced than that seen in either s or Sk. Comparison of values of Sk and Ku* for the final surfaces in each experiment suggests that higher flows generate beds with greater positive skew and kurtosis (Figures 5e (r = 0.73, and p = 0.097) and 4f (r = 0.96, and p < 0.05), where r is the Pearson Product Moment correlation coefficient). Figure 5 Open in figure viewer PowerPoint Variation in standard deviation, skewness, and kurtosis of the bed elevation PDFs for (a–c; key as in Figure 5 a) all the surfaces in each experiment with time and for (d–f) the final surface of each experiment with flow strength. In Figures 5 a– 5 c, the data points are shown for one of the experiments (e1) for reference.

3.4 Variograms The downstream γ(h x , 0) and cross‐stream γ(0, h y ) variograms for both screeded and water‐worked beds approximate the model described above: that is, each variogram can be subdivided into scaling, transition and saturation regions at short, intermediate, and long lags, respectively. By way of illustration, Figure 6 shows the downstream variograms for all the surfaces measured in experiment e2 and for the final surfaces in each experiment. Figure 6 Open in figure viewer PowerPoint h 1x and h 2x and their 95% confidence intervals for the screeded surface. Figure t = 0 and 1200 min for reference and to highlight the initial and final surfaces. Downstream variograms for (a) all the surfaces measured in experiment e2 and (b) the final surfaces in each experiment. In Figure 6 a, the dashed lines indicate estimates ofandand their 95% confidence intervals for the screeded surface. Figure 6 a also shows the data points for= 0 and 1200 min for reference and to highlight the initial and final surfaces. In general, the spline model provided a good fit to the data and generated coherent and consistent estimates of h 1 and H. Considering all experiments, the goodness of fit statistic r2 for the downstream variograms, for example, varies between 0.991 and 0.999 and the 95% confidence limits associated with the estimates of h 1x and H x are of the order 1–2 mm and 6–8%, respectively. Notwithstanding the high r2 values, plots of the residuals for h x ≤ h 1x reveal a convex upward shape, highlighting the presence of nonlinearity in this region of the variogram. This suggests that a scale‐independent scaling region does not perfectly describe the surfaces at short lags and that the parameters h 1 and H need to be interpreted cautiously. Caution also needs to be exercised in respect to the estimates of h 2 since these are highly variable, subject to large errors and difficult to reconcile with the expected roughness characteristics of the surfaces. It might be expected, for example, that the transition region for the screeded beds would be relatively narrow so that h 2 ≈ h 1 . This is because these surfaces should not have exhibited any substantive variations in bed topography larger than the grain scale. However, with the exception of experiment e1 (h 2x = 37 ± 6 mm), estimates of h 2x are many times greater than h 1x and larger, even, than the maximum particle size (72 ≤ h 2x ≤ 342 mm), with 95% confidence intervals of several tens of millimeters (Figure 6a). These high values probably reflect the difficulty in screeding a gravel bed completely flat and illustrate the sensitivity of h 2 to small variations of the medium‐scale topography. In many cases, the sills are also poorly defined, suggesting that the sampling area may be too small. Values of h 2 for the water‐worked surfaces are generally higher, reflecting the development of the mesoscale topography during the experiments (Figure 3). The number of mesoscale forms contained within the test section was, however, low, and estimates of h 2 , therefore, are highly variable. In the light of these difficulties in generating meaningful and generalizable estimates of h 2 , they are not considered further. As expected, the topography of the screeded surfaces was isotropic with a roughness structure and scaling that was consistent with the random arrangement of individual grains (h 1x ≈ h 1y ≈ D s50 ; H x ≈ H y ≈ 0.5). Interestingly, many values of H are slightly in excess of 0.5. This reflects a deviation from a truly random topographic structure and is thought to reflect the presence of some larger surface grains which, due to their greater radii, means that sequences of heights across their tops exhibit positive autocorrelation. Over time, sediment transport forced a series of bed adjustments characterized by increases in bed roughness ((h x , 0)), the scaling region (h 1x ) and the scaling exponent (H x ), as illustrated by the evolution in variogram form (Figure 6a) and by the corresponding temporal variation in h 1x and H x (Figures 7a and 7b). Estimates of (h y , 0), h 1y , and H y also increased during each run. For all water‐worked surfaces (n = 47), little difference is observed between the downstream and cross‐stream estimates of either h 1 ( ) or H ( . The water‐worked surfaces are, therefore, isotropic in h 1 as well as H. A comparison of the variograms for the final surfaces of each experiment does, however, suggest a dependency of surface structure on flow strength (Figure 6b). This dependency is quantified in Figures 7c and 7d which shows that the scaling region (h 1 ) and exponent (H) in both downstream and cross‐stream directions increased from e1 (τ* ≈ 0.024) to e2 (τ* ≈ 0.029), with no further increase at the highest flow strength in e3 (τ* ≈ 0.046). Figure 7 Open in figure viewer PowerPoint h 1 and Hurst exponents H for (a, b; key as in Figure H y > H x , the differences for e1, e1r, and e3r are not statistically significant (p > 0.05). Variation in length scalesand Hurst exponentsfor (a, b; key as in Figure 7 a) all the surfaces of experiment e2 and its replicate e2r with time (downstream direction only) and for (c, d) the final surface of each experiment with flow strength (downstream (solid symbols) and cross‐stream (open symbols) directions). It should be noted that although Figure 7 d implies, the differences for e1, e1r, and e3r are not statistically significant (> 0.05). The evolution of bed surface character is also evident in the 2‐D variograms. As an example, Figure 8 shows the 2‐D variograms for the surfaces of experiment e2 at t = 0, 30, and 1200 min and the variation in the ratio of the downstream and cross‐stream semivariance ((h x , 0)/(0, h y )). The isotropic nature of the initial surface is confirmed by the near‐circular contours in Figure 8a which indicates that the change in γ with lag distance was similar in all directions. As shown in Figure 8d, the ratios γ(h x , 0)/γ(0, h y ) for the initial surface closely approximate unity: the range of correlated data for all initial surfaces of all the experiments (approximately 1 ≤ h ≤ 50 mm; n = 155) is 0.88 ≤ γ(h x , 0)/γ(0, h y ) ≤ 1.14 with 0.95 ≤ γ(h x , 0)/γ(0, h y ) ≤ 1.05 accounting for 70% of all values. The 2‐D variograms for the water‐worked beds (Figures 8b and 8c) also exhibit circular contour patterns indicative of surface isotropy, but only for short lags (approximately h ≤ 4 mm for which 0.88 ≤ γ(h x , 0)/γ(0, h y ) ≤ 1.10 in experiment e2; Figure 8d). At longer lags, the contours are elliptical (Figures 8b–8d), reflecting the development of surface anisotropy [Nikora et al., 1998; Nikora and Walsh, 2004; Mao et al., 2011]. The principal axes of the ellipses are generally aligned in the flow direction (increasing h x ), indicating the development of a flow‐aligned structure. Figure 8 Open in figure viewer PowerPoint t = 0, (b) 30, and (c) 1200 min. (d) Variation in the ratio (h x , 0)/(0, h y ) with lag for all time intervals in experiment e2 (1 ≤ h ≤ 100 mm). Note that the variogram surfaces Figures γ(h x , h y )/s2 for 1 ≤ h x,y ≤ 100 mm (Figures h x ≤ 150 mm, 1 ≤ h y ≤ 100 mm (Figure Two‐dimensional variograms for the surfaces of experiment e2 at (a)= 0, (b) 30, and (c) 1200 min. (d) Variation in the ratio (, 0)/(0,) with lag for all time intervals in experiment e2 (1 ≤≤ 100 mm). Note that the variogram surfaces Figures 8 a– 8 c show the normalized semivariance)/for 1 ≤≤ 100 mm (Figures 8 a and 8 b) and 1 ≤≤ 150 mm, 1 ≤≤ 100 mm (Figure 8 c). Particles < 4 mm comprise, on average, < 7% of the water‐worked, surface grain‐size distributions. The scale of the surface isotropy identified in Figures 8b and 8c (1 ≤ h ≤ 4 mm) approximates, therefore, the smallest grain‐size fraction (4–5.6 mm) that makes up a significant proportion (≈ 8%) of the bed surfaces. Since the point pairs used to calculate γ at this scale are likely to be on the same grain, the small‐scale surface isotropy of the water‐worked beds can be attributed to the surface roughness of individual particles [see also Hodge et al., 2009b]. At slightly longer lags (4 ≤ h ≤ 45 mm), the ellipsoid contours are thought to reflect the preferential orientation of grains with the long axis of the ellipse indicating the alignment of the a axes [Nikora et al., 1998; Nikora and Walsh, 2004; Hodge et al., 2009b]. In this study, surface anisotropy extends beyond the grain scale which suggests that the flow‐aligned structure is some combination of the prevalent alignment of individual particles, clusters of grains and the mesoscale topography described above (Figure 3). Figure 8 confirms that the mesoscale structure developed relatively quickly with the onset of sediment transport and was maintained thereafter for the duration of each run with relatively little modification, despite the changing spatial organization of the component forms (the patches of higher and lower bed elevations).

3.5 Area Ratio As highlighted by Brasington et al. [2012], A R statistics are highly sensitive to DEM resolution, with increasing grid size associated with lower A R values and reduced topographic complexity. In this study, we examine the nature of the adjustment of surface complexity at DEM resolutions of 60 and 5 mm (Figure 9). These scales were chosen on the basis that the results characterize the surface at scales of 180 and 15 mm, respectively, which might, in turn, be expected to correspond with the mesoscale and grain‐scale components of the surface topography (100–200 mm (Figure 3) and 13–17 mm (Table 1), respectively). As anticipated, a comparison of Figures 9a and 9b indicates that the coarser scale surfaces were characterized by less surface complexity (lower Ar values). Figure 9 Open in figure viewer PowerPoint Variation in A R with time for each experiment as measured at DEM resolutions of (a) 60 mm and (b) 5 mm. Since the calculation of A R is undertaken using a 3 × 3 moving window, the surface characterization is at scales of 180 and 15 mm, respectively. Data points are shown for one of the experiments (e1) for reference. At the coarser scale (Figure 9a), the experimental surfaces were initially flat (A R ≈ 0 at t = 0 min) but developed a more complex topography (A R > 0) relatively quickly with the onset of sediment transport. At this scale, the increased complexity of the water‐worked beds is consistent with the development of the flow‐aligned topography during the experiments as seen in the DEMs (compare, for example, Figures 3b, 3e, and 3h). Interestingly, Figure 9a suggests that e2 and e3 developed more complex surfaces than e1, a difference which is also evident in the DEMs (compare Figures 3g–3i). Similar patterns of adjustment in surface complexity were observed at the smaller scale (Figure 9b), and their resemblance to the patterns of grain‐size adjustment shown in Figure 2b suggests that they reflect surface coarsening and the development of a mobile armor.

3.6 Inclination Values of I for the initial and final surfaces of each experiment are compared in Figure 10a. Calculated at the resolution of the DEMs (0.15 mm), I could have been dominated by grain‐surface roughness. However, for the relatively smooth gravels of the experimental mixture, grey‐scale plots with i+ as black and i− as white clearly show grain‐scale patterns, which demonstrate that the index does reflect particle imbrication (i+ and i− defined in Table 2). Figure 10 Open in figure viewer PowerPoint I for (a) the initial and final surfaces in each experiment with flow strength and for (b) all the surfaces in each experiment with time. In Figure Variation infor (a) the initial and final surfaces in each experiment with flow strength and for (b) all the surfaces in each experiment with time. In Figure 10 b, data points are shown for one of the experiments (e1) for reference. Values for the initial surfaces approximate zero (−0.04 ≤ I ≤ 0.04) indicating that the number of positive and negative inclinations was approximately equal, which is the expected result for screeded beds in which particles show no preferential orientation. Interestingly, however, the index is slightly negative for five of the six screeded beds, which indicates a slight bias toward particles with an upstream inclination and may reflect the fact that the beds were screeded in the direction of flow. In contrast, the inclination indices for the water‐worked surfaces vary between −0.14 and −0.09 indicating a predominance of upstream facing slopes and the presence of grain imbrication. The degree of grain imbrication shown by the final surfaces is not sensitive to flow strength (Figure 10a) and is shown to have developed very quickly with the onset of sediment transport (Figure 10b).