Quantifying the Variability in Flow Competence and Streambed Mobility with Water Discharge in a Gravel-Bed Channel: River Esva, NW Spain

Abstract

Streambed mobility in gravel-bed rivers is largely controlled by the rate at which particles with different grain sizes are recruited from the riverbed into the bed load. In this paper, we present a study case in which we explored this question, based on combining field observations using painted plots and the grain size analysis of a large flood sediment deposit in the River Esva, northwest Spain, and the generalized threshold model (GTM) competence model developed by Recking. The main aim was to accomplish a complete characterization of streambed mobility in this river. The obtained results suggest the large potential of the GTM model compared to previous competence models when searching for the quantitative description of particle entrainment and streambed mobility in the River Esva. We observed how the grain size of the bed load in the River Esva tended to be closer to that of the sub-armour bed material during large floods, while moderate magnitude flows tended to carry a relatively fine bed load. Additionally, we compared our results with previously published field observations on flow competence. This comparison outlined the large degree of site specificity in the links between grain size of the bed load and that of the bed material.

1. Introduction

Unravelling the complex interplays between flow strength, grain-size distribution (GSD) of the bed load, and that of the bed material is a central issue in river morphodynamics [1,2]. This topic has been classically approached in fluvial geomorphology and sediment transport studies through the analysis of the thresholds for bed particle entrainment as long as the ability of the stream to displace riverbed particles controls the GSD of the bed load in alluvial channels [2,3,4,5]. In this regard, the inception of motion for streambed particles also exerts significant control over the area of bed surface that is actively involved in bed load motion and the vertical extent of scour and fill processes. Thus, particle entrainment is closely related to the overall streambed stability and the spatial extent of bed disturbance caused by sediment erosion and deposition during floods [6], which ultimately represent important drivers of aquatic habitat quality in river ecosystems [7,8,9,10].

Additionally, particle entrainment is closely related to the “flow competence” concept, introduced first by Gilbert and Murphy (1914) [11] and defined based on the size of the largest clast moved by a given flow. The idea of river competence relies on assuming that the coarsest clast transported by a streamflow provides a measure of the stream’s ability to entrain the bed sediment. Consequently, the measure of the largest clast size of fluvial deposits has been a common and important tool for most geologists. This procedure was frequently used as an instrument to infer velocities, discharges, mean stresses, and depths of flows that transported the particles found in fluvial conglomerates [12,13,14,15,16,17,18,19], or in the evaluation of extreme flood hydraulics [20]. However, the use of the largest moved clasts in this way, i.e., as a technique to infer differences in flow strength and stream’s ability to convey sediment between different sites, involves the implicit assumption that finer grains are always entrained at lower flows than coarser grains and over a relatively broad range of flow discharges [21,22]. This implies ignoring important physical aspects of sediment entrainment that nuance the (in principle) expectable “weight selective” assumption for particle entrainment [2,23,24].

Consequently, there is a clear interest in better understanding the links between the GSD of the bed material, the streambed surface, the bed load, and the spatial–temporal patterns of streambed mobility in coarse bed rivers. Previous works have already explored the relationship between the grain size of the bed load, sampled using sediment traps and/or different samplers, and the grain size of the bed material [1,13,18,25]. However, in many cases, these observations were made before some of the more modern bed load transport and flow competence models were developed (i.e., Parker (1990) [26]; Wilcock and Crowe (2003) [27]; Recking, (2016) [2]), so the former could still benefit from new field testing. Additionally, available field measures of the GSD of the bed load were carried out in many cases during moderate and regular bed load transport episodes (when sampling bed load with Helley-Smith and other kinds of samplers is easier), whereas information about relatively strong events is still scarce. In summary, more field studies are needed in order to better grasp several questions related to streambed mobility, such as understanding whether the GSD of the bed material proxies that of the bed load or not, or how the bed state modulates the relationship between particle entrainment and flow strength.

In this regard, the current paper presents a case study accomplished in a gravel bed stream from northwest Spain, where an extreme flood episode occurred in January 2019. It left important sediment deposits, and its GSD was sampled in the field. Additionally, painted plots were employed in order to monitor streambed mobility during regular, low-magnitude floods occurring before this major flood episode. Both sources of data were used in order to calibrate a streambed mobility model, based on Recking (2016) [2], which, in turn, was used for modelling the GSD of the average annual bed load and the temporal patterns of streambed mobility. The obtained results highlight how the grain size of the bed load in the River Esva tends to be closer to that of the sub-armour bed sediment during large floods, while moderate magnitude flows carry a relatively finer bed load. These results show the potential of coupling a competence model (such as the GTM model [2]) with field observations in order to predict streambed mobility. We also compared our measures and estimates with the bed load information available in Milhous (1973) [28], Andrews (1994) [29], Lisle (1995) [25] and Whittaker and Potts (2007) [1]. This comparison outlined how the relations between grain size of the bed load and that of the bed material can be largely variable with flow strength and among different sites.

2. Rationale

2.1. Physical Basis for Flow Competence Models

Flow competence and streambed mobility are closely related to particle entrainment and the hydraulic thresholds for incipient sediment motion. In this regard, classical research tried to relate shear stress and/or flow velocity to the size of the largest particles set in motion [30,31,32,33]. Initially, threshold stresses for sediment entrainment were related to particle size [30]:

$${\tau }_c={\tau }_c^{*}\times \left({\rho }_s-\rho \right)\times g\times D$$

(1)

where ρ is the density of water, ρ_s that of sediment, g is the gravity acceleration, D the grain size of the bed sediment, τ_c is the critical bed shear stress for incipient motion (dimensional), and τ_c* is the critical dimensionless Shields stress, for which a large range of values have been proposed, normally between 0.03 and 0.06 [23,34]. Once known the value of this parameter, the grain size of the particles that are mobilized with increasing flow discharges (D_max) are easy to estimate:

$$D_{max}=\frac{\tau }{{\tau }_c^{*}\times \left({\rho }_s-\rho \right)\times g}$$

(2)

The previous model Equation (2) holds well in the ideal case of uniform sized particles. In the case of a graded mixture, some works revealed how the effective stress acting on a given grain is strongly influenced by the relative size of the grain within the mixture [35,36]. Sand content seems to influence gravel mobility by reducing the flow resistance and intergranular friction [11,27,37]. Clast arrangements and imbrications may also have a large impact on streambed mobility [38]. As a result, a large variability of situations has been documented regarding the relative mobility of each size class of a sediment mixture in gravel-bed rivers with several studies reporting size selective mobility over a significant range of flow [13,15,28,39,40,41], while others suggested that all sizes in a sediment mixture may begin moving over a range of narrow flow conditions [42,43,44,45,46]. So-called hiding functions were proposed in order to mathematically describe this range of situations:

$${\tau }_{ci}^{*}={\tau }_{cRef}^{*}\times {\left(\frac{D_i}{D_{ref}}\right)}^{\alpha }$$

(3)

where τ_ci is the critical Shields parameter for incipient motion of a given diameter D_i and τ_cref the critical Shields parameter for incipient motion of a reference particle diameter D_ref; α is a coefficient whose value varies from 0 in the case of a total size selective behaviour (threshold conditions only depend on particle size) to −1 in case of perfect equal mobility conditions (all sizes are entrained at the same threshold stress, independent of their size). Rearranging Equation (3), we can estimate again the size of the largest clast moved during floods:

$$D_{max}=D_{ref}\times {\left(\frac{{\tau }_{ci}^{*}}{{\tau }_{cRef}^{*}}\right)}^{\frac{1}{1+\alpha }}$$

(4)

In spite of the great improvement hiding functions (such as Equations (3) and (4)) provide for modelling flow competence and streambed mobility, there are still some open issues. For instance, some researchers have found that a simple power law (Equations (3) and (4)) is not enough to represent the complete relation for any particle size critical stresses in gravel mixtures [34,47]. In addition, gravel-bed rivers are systematically armored [48,49,50,51,52], i.e., the streambed surface is commonly coarser than the underlying sub-armor deposit. The potential changes in the shape of the hiding parameter α with armor break-up are not considered when using a single hiding function. Additionally, the entire bed load GSD cannot be computed directly from Equation (4) (which only supplies the coarsest mobile size) but must be deduced from the bed load rates computed for each size class [26,27,34,49]. This commonly requires the use of an indirect scaling approach [3], implying that the same bed load function can capture (on average) the transport of each individual size class for a given shear stress, an idea that still needs field validation [2].

2.2. Recking’s (2016) GTM Model

Because of some of the issues stated above, Recking (2016) [2] proposed a new flow competence model (named “generalized threshold model”, GTM), as an alternative approach for characterizing streambed mobility and computing the GSD of the bed load. The GTM model computes the GSD of the bed load after modelling which fraction of sediment in each size range is entrained from and which fraction remains static on the riverbed. It does so without assuming a particular shape for the bed load function, in contrast to common approaches based on the hiding function (i.e., Parker (1990) [26]; Wilcock and Crowe (2003) [27]). The GTM workflow operates in two steps. First, it computes the maximum competent size from the following expression (see more details in Recking [2]):

$$M=84\times {\left(\frac{{\tau }_{84}^{*}}{{\tau }_{c84}^{*}}\right)}^{\beta }$$

(5)

where M is the Mth percentile of the surface GSD corresponding to the largest mobile class, τ₈₄* is the dimensionless Shields stress estimated based on the 84th percentile of the surface GSD (D₈₄), τ_c84* is the critical Shields stress for entraining the surface D₈₄. The use of the 84-th percentile as reference parameter in this model is meaningful, insofar as the inception of motion of the coarser particles defines the moment when the armour layer starts to become disorganized [53,54]. Then, a transport stage (τ₈₄*/τ_c84*) close to 1 defines the threshold conditions when the armour layer destabilizes. In this regard, the β exponent establishes the rate at which progressively coarser grains are set in motion with increasing flow strengths. It may take a different value depending on whether or not the armour layer destabilizes. A value of ~0.25 is suggested by Recking [2] when the armour layer is broken (τ₈₄*/τ_c84* > 1) which implies assuming full mobility conditions at τ₈₄*/τ_c84 > 2 as Wilcock and McArdell (1997) [46] reported in the flume. In cases where the armour layer is not broken (τ₈₄*/τ_c84 < 1), there is more uncertainty in choosing the right value for β as it probably depends on the degree of bed armoring/pavement. High β simulates an abrupt increase in competent size when armour breaks (i.e., a well-paved stream), whereas a low value models the case of a bed with a more progressive transition from partial to full mobility conditions.

The second step of the GTM model is calculating the fraction (φ) of bed particles mobilized in individual size ranges (i) of the bed surface:

$${\phi }_i=0.01+\left(1-{\left(\frac{i}{M}\right)}^{\gamma }\right)$$

(6)

where γ can be estimated from:

$$\gamma =0.1+0.5\times {\left(\frac{{\tau }_{84}^{*}}{{\tau }_{c84}^{*}}\right)}^{{\gamma }_2}$$

(7)

where γ₂ is a parameter that needs to be calibrated. The combination of Equations (6) and (7) is very flexible and allows to model a broad range of streambed behaviors, such as progressive mobilizations of surface particles or more abrupt transitions from stability to mobility. Thus, a right γ₂ parameter, based on previous empirical information about riverbed behavior, allows the proper modelling of a given study reach. This is one of the strongest points of the GTM model: the fact of identifying the key parameters (β and γ₂) driving streambed behavior and constraining the competence model. The definition of these parameters from field observations permits to estimate the GSD of the bed load with the following expression:

$$F_i=\frac{{\phi }_i}{{\sum }_{j=i}^M{\phi }_j}$$

(8)

In summary, there are three main advantages of the GTM model compared to previous approaches to flow competence: (1) its flexibility, which allows to model a large range of particle entrainment and streambed stability conditions; (2) the model is well-rooted into our current knowledge on sediment transport mechanics and considers some situations that were not adequately grasped by previous flow competence models; i.e., the change in entrainment conditions with armour break up; and (3) while most of the previous flow competence models focused mainly on computing the mobilized grain sizes, the GTM model also allows the estimation of the mobile fraction in each size range which is at least as important as predicting the maximum transported grain in terms of evaluating streambed disturbance (interesting, for example in many ecological applications). In spite of all these advantages, there have been no previous attempts to take advantage of the GTM model for the characterization of streambed mobility in specific field cases (Vázquez-Tarrío et al. (2019) [55] being the only exception to the best of our knowledge).

2.3. Proposed Workflow

The GTM model provides a well-suited practical framework for streambed mobility studies. Hence, we propose a field workflow adapted to the characterization of streambed mobility in gravel-bed rivers, which combines the flexibility of the GTM model with field observations (Figure 1). It combines two different sources of field data. Painted tracers, which can be employed in order to characterize streambed mobility at low transport conditions, and, as consequence, to obtain the best value for the β exponent in Equation (5) when τ₈₄*/τ_c84* < 1. The second source of data comes from collecting in the field the GSD information of sediment deposited during large floods. This kind of data helps to calibrate Equation (6) by forward modelling: by changing the γ parameter in Equation (6) until obtaining an adequate value to fit the field measures. Finally, the γ parameter obtained from flood deposits could be combined with the field identification of incipient motion conditions in order to calibrate the γ₂ parameter in Equation (7). This paper illustrates the implementation, potential and capabilities of this workflow, using the River Esva (NW Spain) as study case.

3. Study Site

The River Esva is located in the northern watershed of the Cantabrian Mountains, in NW Spain (Figure 2). The climate is temperate-oceanic, with average monthly rainfalls ranging between 52 and 146 mm. The River Esva drains a surface of 464 km². The flood regime is perennial and rainfall dominated, with larger discharges in winter. Average annual discharge is 10.5 m³/s, with average maximum and minimum annual discharges of 19.2 and 4.5 m³/s, respectively.

We selected a ~1 km study reach located close to the village of Trevias (Figure 3). It is a single and straight river reach with an average ~0.003 channel slope and 25 m width. This reach was embanked in 1993 along its entire length in order to protect Trevias from overbank flows.

A lateral gravel bar develops on the left bank in the downstream part of the study reach. Streambed sediment is mainly composed of sub-angular and rod shaped sandstone gravels and cobbles. One square plot of blue-painted stones was deployed on this gravel bar in November 2016 (hereinafter called PP1 painted plot, for more details see Section 4). A second forced bar deposit is located on the left bank in the upstream part of the study reach, where a second plot of yellow-painted stones was deployed in May 2017 (hereinafter named PP2 painted plot, more details in Section 2) (Figure 3).

Additionally, a gauging station located at the downstream part of the study reach (Figure 3) provided water level and flow-discharge records during the whole study period (November 2016 to May 2019) (Figure 4).

4. Materials and Methods

4.1. Painted Plots

As outlined above, two plots of painted stones were deployed in the study site (Figure 3). The painted plots were deployed on the main body of two lateral gravel-bars. As far as our main aim with the plots was monitoring incipient motion conditions for coarse sediment particles, doing so with painted stones is easier in those areas of the channel that are not permanently under water, such as bars.

The first plot of tagged stones (PP1) covers a surface of 1 m × 1 m square on the riverbed. The stones were painted blue in situ, without moving or disturbing in any moment their position. Each grain was labelled and measured with a ruler. Additionally, we took several photographs to record the initial state of the painted stones. In May 2017, 100 stones were collected from the riverbed and taken to the lab, painted yellow, measured, and labelled. We seeded them on the riverbed. Several deployment options were considered for this second generation of tagged stones such as a square plot, a transverse line on the channel bed, or the random distribution of the stones on the streambed. Finally, we decided to seed the stones defining a square of tracers (1 m × 1 m) on the streambed (PP2) (Figure 3). Three were the main reasons: (1) to assure that all the stones were submitted to similar hydraulic conditions; (2) to keep coherence with the strategy previously followed with PP1; and (3) the identification in the field of incipient displacements and stone disturbance is easier with a painted plot

We made regular visits to the field (once every month, on average, and after each major peak of flow), in order to monitor the changes experienced by the painted plots. Each time the painted stones were displaced, we measured the grain size of the largest moved clast. To help this task, we compared photographs of the painted plots taken before and after the transport episode. We also measured the downstream transport distance travelled by the recovered clasts and registered the amount of painted stones recovered and lost after each flood.

4.2. Grain Size Analysis of Flood Deposits

In January 2019, a high-magnitude flow event occurred in the study river. Aerial photographs were available before and after the flood episode, thanks to a monitoring program with an unmanned aerial vehicle (UAV-drone platform). This helped the identification of a conspicuous gravel body deposited on the top of a gravel bar located in the study reach, undoubtedly related to the large flood. This sediment accumulation consists in an elongated gravel-sheet body overlapping a lateral bank deposit along the left margin of the river (Figure 5). Field visits were carried out to identify this deposit and to sample its GSD. We also sampled the GSD of the streambed surface and subsurface.

In the case of the streambed surface and the flood deposit, we used the Wolman (1954) [56] pebble count method. Each pebble count consisted of 200 grains collected along two ~50 m sampling lines spaced ~1–5 m apart. The sampling of grains was done systematically, extracting them at every 1 m intersection along a tape (around twice the largest grain size visually estimated in the field). To minimize the operator’s bias, all the grains were selected and measured by the same person. Metallic templates were used to measure the b-axis of grains > 8 mm and to classify them into half-Ψ size classes. Smaller grains were classified into two groups: grains between 4 and 8 mm and grains < 4 mm.

To sample the GSD of the subarmoured bed material, we used a variant of the pebble count method [56] adapted for the subsurface sediment [57,58]. A surface of 1.5 × 1.5 m was selected on the streambed, and the most superficial and coarse sediment layer was retired. We then mixed the underlying sediment using a shovel. Afterwards, a purpose-designed 10 × 10 cm wide sampling grid was deployed over the 1.5 × 1.5 m square surface. We measured the grain size of each pebble falling in each node of the grid. We repeated this operation at three different emplacements, measuring the size of the b-axis on ~100 subsurface particles per site. Finally, we integrated the three measures into a single GSD for the subarmour bed material. The sampled surface was 100 times larger than the largest clast, so we satisfied the Diplas and Fripp’s (1992) [59] criterion for a representative sample.

4.3. Modelling Streambed Mobility in the River Esva

Data from the painted plots were used with the purpose of defining the maximum mobile grain size at different discharges and to calibrate Equation (5) at low magnitude flows (size-selective conditions). Similarly, the GSD of the January 2019 flood deposit was used for the calibration of Equation (6) at high-magnitude flows (equal mobility conditions). The parameters β (in Equation (5)) and γ (in Equation (6)) were iteratively modified until obtaining those values that minimize the root square mean of differences between the observed (in the field) and the modelled competent sizes (Equation (5)) and GSD of the January 2019 flood deposit (in Equation (6)), respectively. Finally, Equation (7) was calibrated based on the value of γ computed for the January 2019 flood and the γ = 0 assigned to the May 2017 transport episode (close to the incipient motion conditions, see Section 5.3).

The resulting GTM model was applied to characterize streambed mobility in the River Esva. More specifically, we used the calibrated model to answer three questions: i. how the maximum and median competent sizes evolve with flow strength in the Esva river; ii. how the mobile fraction of fine particles (<8 mm) changes with flow stage; and iii. how the number of mobile particles grows with flow discharges. The maximum and median competent sizes were estimated using Equations (5)–(8). The fraction of bed load finer than 8 mm was estimated with Equation (6). The increase in the number of mobile particles was computed based on the mobile fraction of particles in each i-size range:

$$f_s=\frac{\sum {\phi }_i\times p_i}{\sum p_i}$$

(9)

where f_s is the portion of the streambed particles that are entrained and p_i the relative fraction of each i-size class in the surface GSD.

Applying the GTM model requires the computation of shear stresses from flow discharge. To achieve this, we used the Rickenmann and Recking (2011) [60] fit to the Ferguson (2007) [61] flow resistance equation and estimated flow depths (d) from flow discharge using:

$$d=0.015\times D_{84}\times \frac{Q^{* 2p}}{p^{2.5}}$$

(10)

where:

$$Q^{*}=Q/w\times \sqrt{gS{D}_{84}^3}$$

(11)

where Q is the water discharge, Q* the dimensionless water discharge, w the channel width, S the channel slope, and p = 0.24 if Q* < 100 and 0.31 otherwise. Then, assuming a rectangular cross-section, we estimated the hydraulic radius (R) from:

$$R=\raisebox{1ex}{$d\times w$}\!\left/ \!\raisebox{-1ex}{$\left(2d+w\right)$}\right.$$

(12)

and, finally, shear stresses were directly computed from the hydraulic radius–slope product:

$$\tau =\rho \times g\times S\times R$$

(13)

The GTM model also needs a value for the critical threshold of motion of the D₈₄ (τ_c84*). We based it on the slope dependent relation proposed by Recking (2009) [62]:

$${\tau }_{c_{84}}^{*}=\left(1.32\times S+0.037\right)\times {\left(\frac{D_{84}}{D_{50}}\right)}^{-0.93}$$

(14)

Finally, we computed the average GSD of the bed load mobilized during the entire study period (November 2016–July 2019). We used the 15 min discharge record available from the gauging station for the study period. Then, we applied the calibrated GTM model and Equation (8) to estimate the evolution of the GSD of the bed load at each 15 min time step of the flood hydrograph.

Then, the time integrated bed load volumes (V_i) in each size class, during the entire study period, can be estimated from:

$$V_i={\displaystyle \sum }_{T=T_0}^{T=T_f}q{s}_{T·}p{p}_{i_T}\times T$$

(15)

where T₀ and T_f are the departure and final times of the study period (respectively), qs_T is the bed load rate at each time step T, and pp_iT is the proportion of the i-size class in the bed load. Bed load rates were estimated using Recking’s (2013) [63] sediment transport equation, which appears to work well in coarse-bed rivers with the same range of slope and bed conditions as the River Esva [63,64]. Once we computed V_i for the different size classes, we calculated the percentage of each size class in the total bed load volumes (V_T) mobilized during the study period:

$${\%}_i=\raisebox{1ex}{$V_i$}\!\left/ \!\raisebox{-1ex}{$V_T$}\right.$$

(16)

4.4. Comparison to Previous Studies

We compared the GSD of the flood deposit sampled in the field to both the surface and subsurface GSD, in order to see how the grain size of the carried load relates to the grain size of the bed material. Later, in order to nourish the discussion on how the grain size of the carried load may or may not correlate to flow strength and the GSD of the streambed material, we re-analyzed previously published information on the grain size of the bed load available in Milhous (1973) [28], Andrews (1994) [29], Lisle (1995) [25] and Whittaker and Potts (2007) [1] (

Table 1. Data of grain size of the bed load, compiled from previous studies and compared to our own data for the River Esva. S: Slope. W: Channel width (in m). D_50s and D_50ss: Median size of the surface and subsurface grain-size distribution (in mm), respectively. D_50bl: Median size of the bed load. D_max: Largest mobile grain (in mm). Q_bkf: Water discharge at bankfull (m³/s). Q: Water discharge (m³/s) during the sampled sediment transport episodes. Sample: Bed load sediment transport method.

River	Source	S	W	D50s	D50ss	D50bl	Dmax	Qbkf	Q	Sample
Oak creek	[28]	0.0094	3.7	54	20	3–27	8–98	3.4	0.2–3.4	Vortex sampler
Sagehen creek	[29]	0.0102	2.6	58	30	12–65	31–126	0.3	0.2–0.6	HS
Bambi creek	[25]	0.0082	3.6	-	9	3	22	1.7	2.5	HS
Goodwin creek	[25]	0.0033	12.9	12	8	8	24	90.0	81	HS
Jacoby creek	[25]	0.0063	17.2	27	8	4	25	9.0	34.2	HS
NF Caspar creek	[25]	0.0130	4.4	15	9	4	80	3.1	6.2	Pit
Redwood creek 1	[25]	0.0140	44.6	5	4	1	30	430.0	374.1	HS
Tanana river	[25]	0.0008	315.0	13	1	17	22	-	2040	HS
Tom Mac Donald Creek	[25]	0.0060	6.1	19.8	11	8	22	3.6	12.6	HS
Turkey brook	[25]	0.0086	3.0	22	16	11	100	13.0	16.9	Pit
Dupuyer creek	[1]	0.0100	8.0	54	-	47–60	88–155	6.5	4.6–10.3	Basket sampler
River Esva	This study	0.0030	25.2	59	22	-	16–208	-	31.3–183.9	Painted stones

) and we compared them to our field observations.

5. Results

5.1. Painted Plots

Four floods were monitored with the painted stones from November 2016 to March 2019 (

Table 2. Summary of results of the field surveys in the River Esva accomplished here.

Date	Painted Plot	Q (m3/s)	τ*/τc*	Dmax	% Recovery	L (m)
17/01/2017	PP1	41.93	0.73	54.5	86	~0.2
02/03/2017	PP1	74.63	0.98	90	82	4.0
29/05/2017	PP1	31.34	0.63	0	82	0.0
12/11/2017	PP1	183.95	1.55	208	82	35.2
29/05/2017	PP2	31.34	0.63	16	96	<0.1
12/11/2017	PP2	183.95	1.55	158	3	19

Q: Peak flow discharge. τ*/τ_c*: Transport stage ratio (at peak discharge). D_max: Maximum mobile size (mm). L: Maximum distance travelled by the retrieved painted stones.

).

The first of the studied flood episodes occurred on 17 January 2017. It had a 41.93 m³/s peak discharge which corresponds to a frequent, regular peak of flow (<1 year return period discharge). This episode was able to move 23% of the seeded painted stones, but displacements were not very important, i.e., the maximum measured transport distances were relatively low (~0.2 m). The b-axis of the largest moved clast was 55 mm. These field observations can be considered as reliable, as approximately 85% of the tracers initially painted were retrieved. They reflect a partial mobility regime and close to the threshold of motion for the median-sized particles in the streambed.

The second transport episode, which occurred on 2 March 2017, was also monitored with the PP1 painted plot. It corresponded to a 74.63 m³/s discharge (<1 year return period flow). During this peak of discharge, a larger amount of painted tracers was disturbed (55%), but the displacements were again not very important; the maximum measured downstream transport distance was around 4 m. The size of the largest moved clast was 90 mm, coarser than in the previous episode. All these observations suggest partial mobility conditions and bed shear stresses close to the threshold of motion for the coarser particles in the riverbed. We recovered a large number of the initially painted stones (~82%), so all our observations were reliable.

The third transport episode took place on 29 May 2017. It was a 31.34 m³/s peak of discharge that was monitored with both the PP1 and PP2 painted plots. The PP1 painted stones did not experience any apparent movement. Conversely, the stones from the PP2 painted plot experienced some kind of reorganization, consisting in rolling and pivoting of some clasts and a very short downstream dispersion (<0.1 m). The PP2 was painted only six days before this event happened, while the PP1 painted stones had already experienced two previous peaks of flow. Consequently, differences in clast structuration may explain why the positions of PP2 stones were reorganized while those of PP1 clasts remained undisturbed.

The fourth and relative high-magnitude transport episode occurred on 12 November 2017 and was surveyed with both PP1 and PP2 painted stones. It had a peak discharge of 183.95 m³/s, which corresponds to a ~2 year return period flow. This episode displaced 86% of the painted stones from PP1, with some of the retrieved tracers experiencing downstream displacements as far as 35 m. The b-axis of the largest moved clast was 210 mm. We were able to recover 82% of the initially PP1 painted stones, so again the observations made with PP1 tracers could be considered as trustworthy. In contrast, we could only recover three of the initially painted stones from PP2; the rest were lost. The largest clast of the three PP2 retrieved stones measured 158 mm.

5.2. GSD of a Large Flood Deposit

An exceptionally large flood episode occurred in the River Esva between the 19th and 29th of January 2019. Peak discharge during this event reached 545 m³/s, which according to the available gauging records corresponds to a ~50 year return period flow discharge. A large sediment deposit was accumulated during this flood over a gravel bar in our study site (Figure 5). This deposit was first identified with the help of a UAV drone. Later, we recognized it in the field and we measured its GSD using the Wolman (1954) [56] method (Figure 6).

We also compared the GSD of the flood deposit to those of the surface and subsurface streambed material (Figure 6 and

Table 3. Comparison between the grain-size distribution (GSD) of the surface sediment, the subarmour bed material, and the January 2019 flood deposit in the River Esva. Values between brackets correspond to the grain size metrics estimated for the truncated (<8 mm) GSD.

GSD	% Particles < 8 mm	D16 (mm)	D50 (mm)	D84 (mm)
Surface	0	34 (34)	59 (59)	107 (107)
Subsurface	8	11 (20)	22 (30)	47 (58)
January 2019 (flood deposit)	22	3 (20)	22 (30)	42 (49)

). The GSD of the flood deposit was very similar to that of the subsurface GSD, and both were considerably finer than the surface GSD. The similarity between the GSD of the January 2019 flood and the subsurface GSD suggests that the GSD of the bed material was mainly sculpted by large floods along the River Esva. Nevertheless, the flood deposit showed a larger percentage of sediment finer than 16 mm than the subsurface GSD. This difference might be related to winnowing of fine sediment during low floods, coarsening the topmost layer of the streambed and washing relatively fine particles out of the riverbed. However, potential bias when sampling GSD of the subsurface bed material cannot be discarded. Actually, when the three GSD are truncated at 8 mm, the GSD of the January 2019 flood deposit and the subsurface bed material seem to collapse into a very similar curve. This reinforces the idea that the subsurface sediment is depleted in fine sediment compared to the bed load carried during large floods.

5.3. Calibrating the GTM Model

The results from our field observations were used to calibrate the GTM competence model [2] for the River Esva. Painted plots allowed the calibration of Equation (5) which describes how the maximum mobile percentiles evolve with flow stage. We obtained an optimum ~3 − β parameter for transport stages below 1. In the case of transport stages larger than 1, we noticed how the painted plot data were in accordance with the 0.25 − β parameter initially suggested by Recking (2016) [2] that involves full mobility conditions when transport stages are larger than 2 (Figure 7A).

To calibrate Equation (6) (quantifying how mobile fraction in each size evolves with flow strength) at high flows (equal mobility conditions) we based on the GSD of the January 2019 flood deposit. We searched for the value of the γ parameter (in Equation (6)) that best minimizes the root mean square of differences between observed and modelled GSD (Figure 7B). We obtained a 2.58 − γ value. Nevertheless, according to Reference [2], the value of the γ parameter may vary with flow strength (Equation (7)). In order to better constrain our calibrated model, based on the data of 29 May peak of flow (when painted tracers practically did not move), we considered that the fraction of mobile particles in all size classes decreased to almost 0 at transport stages close to 0.6; this involves a γ parameter ~0. According to this and the ~2.6 γ value for transport stages close to 3.3, the best fit for Equation (7) was obtained with a 1.60 γ₂ parameter (Figure 7C).

5.4. Characterizing Streambed Mobility in the River Esva

Using the GTM model calibrated for our study site, we accomplished an analysis on how streambed mobility evolves in the River Esva with flow stage. We estimated four parameters: (1) the median size of the entrained particles; (2) the size of the coarsest mobile clast; (3) the fraction of fine sediment (<8 mm) in the bed load at every flow stage; and (4) the number of bed surface particles that are entrained and actively incorporated into bed load motion with increasing discharges.

The modelled values clearly show how all these four parameters changed with flow strength (Figure 8). There was a remarkable change around ~75 m³/s where both the grain size and the number of entrained particles experienced a sudden jump. This discharge corresponded to a transport stage close to 1; that is to say, close to the moment when the coarser clasts in the riverbed started to incorporate into motion and the armour layer destabilized and broke up.

5.5. Modelling the Long-Term Averaged GSD of the Bed Load

Using the flow discharge records available for the study site, we estimated the GSD of the bed load integrated for the whole study period. This covers almost three entire hydrological years (2016–2017, 2017–2018, and 2018–2019), so we could consider it somewhat a proxy of the annual average GSD of the bed load. The obtained GSD is considerably finer than the GSD of the January 2019 flood deposit and the subsurface bed material (Figure 6). Indeed, this difference persists even when the GSDs are truncated at 8 mm (Figure 6B), so this is not just a consequence of a larger amount of fine sediment: the time integrated GSD of the bed load is still finer than the subsurface bed material, even when the fine sediment bed material load is not considered.

5.6. Comparison to Previous Studies

In Figure 9 and Figure 10 we plotted the River Esva data together with data from previous studies [1,25,28,29]. This comparison shows a statistically significant (p-value < 0.05) correlation between the largest clast size of the bed load and the transport stage ratio. The grain size of the largest transported clast approximates the diameter of the D₈₄ when transport stages were close to 1 (Figure 9A and Figure 10A). This confirms how once the thresholds of motion of the coarser grains are reached, these coarse particles become incorporated into the bed load mass.

Conversely, in the case of the median sizes of the bed load, we did not observe any connection with the flow strength (Figure 9B and Figure 10B). Data from the different study sites appear well segregated which outlines the high degree of site specificity in the relations between GSD of the bed load and flow strength. This highlights the importance of site-specific controls on streambed mobility. Additionally, data plots corresponding to transport stages higher than 1 show an important amount of scatter. This brings into question the idea of a simple relation between the GSD of the subarmour material and that of the bed load during large floods.

Moreover, correlation between the size of the largest clast and flow strength is stronger when the largest moved clast sizes are normalized by the D₈₄ of the surface sediment than when they are normalized using the D₈₄ of the subsurface bed material. This may suggest how the remobilization of the coarser grains from the riverbed surface may largely control the size of the largest clast in motion.

6. Discussion

6.1. General Discussion

Our work illustrates the great potential and interest of combining a simple and open model such as the GTM [2] with relatively easy-to-take field measures in order to characterize the largest mobile size, fine sediment mobility or overall streambed mobility/stability behavior (Figure 8). Indeed, our field observations were in good accordance with the general shape of the GTM model [2] which is partially due to the large flexibility of the GTM model compared to previous competence models (Figure 11).

One main advantage of the GTM model relies is in its ability to model armour break-up and its adaptability to a broad range of armour destabilization conditions. The Esva calibrated GTM model shows two break points at ~75 and ~275 m³/s (Figure 11), corresponding to the flow discharges when the coarser grains started to be entrained and when full mobility conditions were reached, respectively. According to Figure 11, the GTM was able to model a situation where we have size-selective entrainment at low flows, a progressive armour destabilization at moderate flows and a sudden final breakup of the armour layer at very large flows. In fact, this involves some complex combination of size-selective and equal mobility conditions, very suitable for the River Esva (Figure 11). Previous competence models based on power-law hiding functions (Equation (3)) seem unable to grasp such behavior (Figure 11), as they are only able to model a progressive mobilization of the different grain sizes or a sudden mobilization of all the grain sizes at a narrow range of flows. The Komar (1987) [13] fit to Carling (1983) [15] data (involving a α-parameter ~0.85) was the one closest to the calibrated GTM model.

In this work, we did not consider the possibility of temporal changes in surface GSD when calibrating of the GTM model. However, changes in surface GSD could be expected in gravel-bed rivers, driven by variable flow discharges and temporal fluctuations in sediment supply [51,65]. Nevertheless, we accomplished several pebble counts (~100 counts per sample) in the study reach all along the time period of study (Figure 12). Between 2016 and the January 2019 flood, we did not observe significant grain size changes: we observed a ~15% variability around the mean value for the D₅₀ which is comparable to the precision expected for pebble counts based on 100–200 clast-size measures [58,66]. Thus, we could consider that surface GSD remained quite stable between 2016 and 2019. However, after the January 2019 flood, we reported significant streambed fining. This was probably related to the liberation of fines after widespread armour disorganization during this high-magnitude flow.

Consequently, major floods in the River Esva seem to be associated with a temporal replenishment in fines of the streambed surface. Hence, we could expect differences in the GSD of the bed load after major floods, which will persist until the streambed exhausts again the fine sediment. Thus, we repeated the GTM model based on the surface GSD measured after the January 2019 flood and we compared to our previous estimations (Figure 13). In general, we observed how both the median grain size and the largest mobile clast during moderate flow discharges are finer after the January 2019 flood than before; conversely, differences tend to disappear for higher water discharges. On this point, seasonal and inter-annual changes in bed load rating curves have been documented by several authors in gravel-bed rivers [67,68]. We believe that the GTM model has the potential to grasp this kind of seasonal and annual bed load hysteresis patterns. Changes in in-channel sediment supply conditions should have an imprint on bed surface texture, grain size and streambed structuration [38,66]. Hence, if these changes in surface GSD are adequately monitored in the field, then they could be accounted for in the workflow here proposed (Figure 1) and the hysteresis behavior in flow competence could be approached based on the GTM model.

6.2. GSD of the Bed Load and the Streambed Sediment: Implications for Paleo-Hydrological Analysis

According to our results, the GSD of the time-integrated bed load is considerably finer than the GSD of the subsurface bed material. In addition, the comparison between the GSD of the January flood deposit and the GSD of the subsurface bed material suggests that the GSD of the subarmour material proxies that of the bed load carried during large floods in the River Esva. Church and Hassan (2005) [69] also compared the grain size of the bed load to that of the subsurface bed material in Harris Creek (Canada), based on bed load volumes collected with a pit trap. They also observed how the grain size of the bed load approaches that of the subsurface bed material in periods of large floods, and a time averaged bed load finer than the subsurface sediment.

Low-magnitude peaks of flows and recession limbs of major flows seem able to mobilize a large amount of fine sediment, which passes over the Esva’s riverbed and is not stocked into the bed material. Additionally, during large floods, much of the sand fraction of the bed material can escape as suspension load. Both issues help explain why the GSD of the time integrated bed load is finer than the GSD of the sub-armour sediment, although some potential methodological bias cannot be neglected. We defined the GSD of the time integrated bed load based on bed load rates computed using Recking’s (2013) [63] (updated by Recking et al. (2016) [70]) bed load equation, which is said to perform well in gravel-bed rivers [64]. However, at low flows, fine sediment exhaustion from the riverbed surface can give rise to a supply-limited situation, compromising the applicability of bed load formulae. Thus, the amount of fine sediment in the time integrated bed load may be overestimated.

In summary, the GSD of the bed load during large floods tends to resemble the GSD of the subarmour bed material. However, re-analysis of previously published data suggests that this is not a general trend in gravel-bed rivers, and it could be largely site specific. A wide scatter is present in the available data (Figure 9 and Figure 10) and the different data sets appear well-segregated, so site-specific conditions control the existing relations between the GSD of the subsurface bed material and that of bed load during large floods. Consequently, patterns of armour development and destabilization are largely site-specific, influencing the links between GSD of the bed load and that of the bed material.

In this regard, the size of the largest clast has been a classical tool in the evaluation of flow competence by many geologists and sedimentologists interested in the paleohydrological and paleohydraulical interpretation of old fluvial conglomerate deposits, who used it (for instance) for paleo-velocity or paleo-slope analysis [12,13,14,15,16,17,18,19,20,71,72,73,74,75]. However, translating flow competence measures, based on the largest clast size, into paleohydraulic conditions is not direct; as we have stated above, there is no single and general function linking the grain size of the bed load to flow strength. In addition, particle weight is not the only control on particle mobility, but packing, imbrication, clast interlocking, fine sediment content and a broad range of exposure and hiding effects distort the links between particle size and entrainment conditions. On this point, some authors documented equal mobility conditions for all the grain sizes present in the bed, whilst others have reported a much more size-selective behavior.

Moreover, gravel-bed rivers are systematically armored [48,49,50,51,52], i.e., the streambed surface is commonly coarser than the underlying subarmour deposit. In this regard, the topmost layer of the streambed represents the local source of sediment particles during bed load transport and surface coarsening results from the adjustment of the streambed’s surface to the bed load transfers during the dominant channel-forming flows. Thus, the differences between the surface and subsurface GSDs may potentially provide some interesting clues for evaluating streambed mobility and/or paleo-hydrological analysis. Unfortunately, the topmost layer of the streambed is rarely preserved in old fluvial deposits [75]. Consequently, the grain size of the subsurface sediment is the only available information for the stratigraphic application of flow competence. Additionally, grain size of fluvial deposits is probably driven mainly by hydraulic conditions during deposition [76,77] rather than by particle entrainment.

With the above in mind, some doubts arise about the stratigraphic interpretation of flow competences based on the largest grain measured in subarmour deposits, whether they may be a good proxy of large magnitude and less frequent transport episodes, or instead they provide an average picture of the more frequent channel-forming flows. Subsurface grain-size distribution has been commonly considered a good proxy for the grain size of the average annual bed load [53,78,79], but our results suggest that there is not a straightforward relation between the GSD of the subsurface material, the time integrated GSD of the bed load and/or that of the bed load during moderate to large floods. These relations are largely controlled by the specific streambed conditions (armoring/pavement, packing, imbrications, etc.) at a given river reach, so the paleo-hydrological interpretation of large clast measures in fluvial conglomerates could not be straightforward.

6.3. Frequency and Intensity of Streambed Mobility: Implications for Streambed Stability Analysis

The timing and intensity of riverbed mobility are of great interest in river ecology [6], since they determine the evolution of bed surface structure, which is the physical support of habitat for many aquatic organisms. Therefore, it is not surprising that stream ecologists have also used the flow competence concept in order to evaluate stream stability during flows [8,22,80,81] in relation to several issues, such as streambed perturbation [6] or the response of benthic macroinvertebrates to floods [82,83,84,85].

The magnitude and frequency of coarse particle motion contributes to shape habitat for macroinvertebrates that live on the riverbed surface and defines the bed conditions that render river gravel more stable through time between large floods. A reductionist view of gravel riverbed mobility based on the flow competence concept would suggest that the timing and intensity of streambed disturbance is driven by the ability of the channel to set in motion the coarser clasts. However, this overlooks the importance of other questions. For instance, the fact that a clast in a given class size is entrained does not mean that all clasts in that size are actually moving (partial mobility conditions [86]). We believe that the workflow outlined here is well suited for the characterization of streambed mobility in gravel-bed rivers, always bearing in mind the site-specific character of each riverbed.

Our results show an abrupt change in the streambed behavior at discharges able to destabilize the armored surface. The size of the largest moved clast and the median size of the bed load show a sudden increase. Furthermore, the percentage of fine sediment in the bed load decreases. This relates to a larger amount of coarser clasts set in motion that decreases the relative proportion of fine sediment. The percentage of mobile particles in the riverbed also increases, reaching a limiting value of ~80% during large floods as full mobility conditions are entirely never reached. This implies that an important number of clasts remain static during floods. This contrasts with observations made by Lisle et al. (2000) [87] and Haschenburger and Wilcock (2003) [88] who found that only a small portion of the bed is inactive during bankfull floods. Nevertheless, other authors have also observed that a large amount of the riverbed could be immobile during large floods [6,89,90], in agreement with our observations.

7. Conclusions

In this work, we combined field observations of particle entrainment using painted plots and grain size analysis of a large flood deposit, with the GTM model, recently developed [2], in order to quantify streambed mobility in the River Esva. The GTM model proved to be a very flexible tool for the characterization of patterns of particle entrainment and armour destabilization in our study site, and we illustrated the potential of this workflow to address these complex issues.

Additionally, we compared our results with field observations from previous studies that we compiled from the literature. We observed how the grain size of the bed load tended to be closer to that of the sub-armour bed material during large floods, while moderate magnitude flows tended to carry a relatively fine bed load. However, there was a huge scatter in the data. Hence, inferring the grain size of the bed load or paleo-flow conditions from the grain size analysis of flood deposits should be done with care.