Emeline Perret
email: e.perret@cnr.tm.fr

Benoit Camenen
email: benoit.camenen@inrae.fr

Céline Berni
email: celine.berni@inrae.fr

Kamal El Kadi Abderrezzak
email: kamal.el-kadi-abderrezzak@edf.fr

Benjamin Renard
email: benjamin.renard@inrae.fr

Uncertainties in models predicting critical bed shear stress of cohesionless particles

Our data show a large scatter for the critical Shields stress for initial sediment motion. The main sources of dispersion are related to the methodological procedures defining the inception of movement (i.e., visual observations or extrapolation of sediment transport rate ) and to the estimation of the bed shear stress. The threshold for sediment motion varies with many factors related not only to grain size, but also with bed composition (e.g., presence of fine sediments in a coarse matrix), arrangement (e.g., bed roughness, grain orientation and characteristic lengths 1

Perret, September 2, 2022 of bed structures) and slope. New models to estimate the critical Shields number are proposed combining both grain size or/and bed slope. Model parameters and uncertainty are estimated through Bayesian inference using prior knowledge on those parameters and measured data. Apart from the uncertainty in observations, two types of uncertainty can be evaluated: one related to the parameter estimation (i.e., parametric) and one related to the choice of the model (i.e., structural).

Eventually, a four-parameter model based on both the grain size and bed slope yields the best results and demonstrates a potential interaction between these two parameters. Model uncertainty remains, however, large, which indicates that other input parameters may be needed to improve the proposed model.

INTRODUCTION

Understanding sediment transport is a major concern in many fluvial and ecohydraulic studies (e.g., riverbed mobility, habitat, water quality) and predicting the critical conditions for incipient particle motion remains a fundamental and practical problem. Bedload increases rapidly and nonlinearly with bed shear stress, and large uncertainties in predicting its rate near incipient motion have been observed in gravel-bed rivers [START_REF] Camenen | A bedload sediment transport formula for the nearshore[END_REF][START_REF] Recking | Feedback between bed load transport and flow resistance in gravel and cobble bed rivers[END_REF][START_REF] Camenen | Assessment of methods used in 1D models for computing bed-load transport in a large river: the Danube River in Slovakia[END_REF]. [START_REF] Shields | Anwendung der Ahnlichkeits-mechanik und der turbulenzforschung auf die geshiebebewegung [application of similarity principles and turbulence research to bed-load movement[END_REF] defined the dimensionless bed shear stress as:

𝜃 = 𝜏 (𝜌 𝑠 -𝜌)𝑔𝑑 (1)
with 𝜏 the bed shear stress, 𝜌 𝑠 and 𝜌 the densities of sediment and water, respectively, 𝑔 the acceleration of gravity, and 𝑑 the grain size. The criterion for incipient motion of sediment particles is commonly expressed in terms of the critical Shields number 𝜃 𝑐𝑟 . Most sediment transport formulas, generally derived from laboratory experiments on well-sorted sediment mixtures, relate bedload rate 𝑞 𝑠𝑏 to the excess bed shear stress (𝜃 -𝜃 𝑐𝑟 ) [START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF][START_REF] Parker | Bed load and size distribution in paved gravel-bed streams[END_REF][START_REF] Van Rijn | Sediment transport, part I : bed load transport[END_REF][START_REF] Lajeunesse | Bed load transport in turbulent flow at the grain scale: Experiments and modeling[END_REF]. The validity of these formulas may be questionable when applied to field cases, such as gravel bed rivers with poorly sorted sediment mixtures and complex bed features [START_REF] Recking | A comparison between flume and field bedload transport data and consequences for surface based bedload transport prediction[END_REF]. Accurate estimation of the bed shear stress and its critical 2 Perret, September 2, 2022 value for incipient motion is then challenging [START_REF] Perret | How does bed surface impact bedload transport over gravel-bed rivers?[END_REF]. [START_REF] Buffington | A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers[END_REF] reported a large dataset for the critical Shields number and the [START_REF] Soulsby | Threshold of sediment motion in coastal environment[END_REF] equation provides a rough fit of 𝜃 𝑐𝑟 expressed in terms of the dimensionless grain size 𝑑 * = 𝑑 50 [𝑔(𝑠 -1)/𝜈 2 ] 1/3 (with 𝑠 = 𝜌 𝑠 /𝜌 the relative sediment density, 𝜈 the kinematic viscosity of the fluid, and 𝑑 50 the median grain size). Still, a significant scatter in the data exists, as for a given 𝑑 * -value, 𝜃 𝑐𝑟 can vary more than one order of magnitude. Data scatter may result from the experimental set-up conditions (e.g., initial bed arrangement) and from the methodological procedures used to define the concept of incipient motion and to compute 𝜃 𝑐𝑟 [START_REF] Buffington | A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers[END_REF]. The scatter in the data may also reflect that 𝜃 𝑐𝑟 depends not only on grain size [START_REF] Garcia | Sedimentation Engineering: Processes, Measurements, Modeling, and Practice[END_REF]) but also on bed slope [START_REF] Recking | Theoretical development on the effects of changing flow hydraulics on incipient bedload motion[END_REF], hiding/exposure of grains [START_REF] Wilcock | Surface-based transport model for mixed-size sediment[END_REF], particle imbrication, and degree of clogging [START_REF] Perret | Transport of moderately sorted gravel at low bed shear stresses: the role of fine sediment infiltration[END_REF].

Several studies have put forward the dependence of the critical bed shear stress on bed arrangement [START_REF] Tait | The physical processes of bed armouring in mixed grain sediment transport[END_REF][START_REF] Haynes | Stress history effects on gravel bed stability[END_REF][START_REF] Yager | Resistance is not futile: Grain resistance controls on observed critical shields stress variations[END_REF][START_REF] Perret | How does bed surface impact bedload transport over gravel-bed rivers?[END_REF][START_REF] Hassan | Experimental insights into the threshold of motion in alluvial channels: Sediment supply and streambed state[END_REF][START_REF] Hodge | X-ray computed tomography reveals that grain protrusion controls critical shear stress for entrainment of fluvial gravels[END_REF]) which has been described through many indicators, such as the roughness height of grains, their shape [START_REF] Lane | Some observations on the effect of particle shape on the movement of coarse sediments[END_REF][START_REF] Li | Laboratory measurements of pivoting angle for applications to selective entrainment of gravel in a current[END_REF][START_REF] Petit | Evaluation des critères de mise en mouvement et de transport de la charge de fond en milieu naturel[END_REF], emergence [START_REF] Fenton | Initial movement of grains on a stream bed: the effect of relative protrusion[END_REF], orientation and imbrication [START_REF] Laronne | Interrelationships between bed morphology and bedmaterial transport for a small, gravel-bed channel[END_REF][START_REF] Reid | The continuous measurement of bedload discharge[END_REF][START_REF] Brayshaw | The hydrodynamics of particle clusters and sediment entrainment in coarse alluvial channels[END_REF], the degree of bed armouring, and the characteristic lengths of bed clusters/ structures [START_REF] Church | Stabilizing self-organized structures in gravelbed stream channels: Field and experimental observations[END_REF][START_REF] Venditti | Bed load bias: Comparison of measurements obtained using two (76 and 152 mm) Helley-Smith samplers in a gravel bed river[END_REF]. Because the antecedent flow conditions impact the arrangement of the bed surface, 𝜃 𝑐𝑟 is thus related to the stress history [START_REF] Haynes | Stress history effects on gravel bed stability[END_REF]. The critical Shields number of coarse particles can also vary by several percent [START_REF] Perret | Transport of moderately sorted gravel at low bed shear stresses: the role of fine sediment infiltration[END_REF] according to the proportion of matrix fines (cohesive or not) [START_REF] Reid | The incidence and nature of bedload transport during flood flows in coarse-grained alluvial channels[END_REF][START_REF] Curran | The decrease in shear stress and increase in transport rates subsequent to an increase in sand supply to a gravel-bed channel[END_REF][START_REF] Jain | Cohesion influences on erosion and bed load transport[END_REF][START_REF] Barzilai | Effect of changes in fine-grained matrix on bedload sediment transport in a gravel-bed river[END_REF][START_REF] Kuhnle | Sand transport over an immobile gravel substrate[END_REF][START_REF] Wren | Turbulent flow and sand transport over a cobble bed in a laboratory flume[END_REF][START_REF] Perret | Transport of moderately sorted gravel at low bed shear stresses: the role of fine sediment infiltration[END_REF], i.e., fine sand can increase bedload by lubrication, whereas the opposite effect is observed with silt and clay due to consolidation effect. Finally, hiding and exposure modify the critical Shields number for each size class in mixtures of non-cohesive sediment particles [START_REF] Jackson | Influences of increased sand delivery on the morphology of sand and gravel channel[END_REF][START_REF] Ikeda | Experimental study of heterogeneous sediment transport[END_REF][START_REF] Wilcock | Surface-based transport model for mixed-size sediment[END_REF][START_REF] Curran | The decrease in shear stress and increase in transport rates subsequent to an increase in sand supply to a gravel-bed channel[END_REF][START_REF] Kuhnle | Sand transport over an immobile gravel substrate[END_REF][START_REF] Wren | Turbulent flow and sand transport over a cobble bed in a laboratory flume[END_REF][START_REF] Perret | Transport of moderately sorted gravel at low bed shear stresses: the role of fine sediment infiltration[END_REF]). Nevertheless, the effects of hiding/exposure can be quantified 3 Perret, September 2, 2022 based only on a reference critical bed shear stress for unisized material or based on the median grain size -the focus of the present study.

Several have observed that 𝜃 𝑐𝑟 increases with mild slopes longitudinal bed slope 𝑆 (0.001 ≤ 𝑆 ≤ 0.05) [START_REF] Shvidchenko | Flume study of the effect of relative depth on the incipient motion of coarse uniform sediments[END_REF][START_REF] Mueller | Variation in the reference shields stress for bed load transport in gravel-bed streams and rivers[END_REF][START_REF] Lamb | Is the critical shields stress for incipient sediment motion dependent on channel-bed slope?[END_REF][START_REF] Recking | Theoretical development on the effects of changing flow hydraulics on incipient bedload motion[END_REF].

For very steep slopes (𝑆 > 0.05), [START_REF] Chiew | Incipient sediment motion on non-horizontal slopes[END_REF] demonstrated that 𝜃 𝑐𝑟 decreases with 𝑆.

The reasons for the increase in 𝜃 𝑐𝑟 with 𝑆 ≤ 0.05 remains partially explored. When 𝑆 increases, stable bed structures appear, leading to morphologic changes and less available shear stresses for bedload. The slope effect could in fact be a drag effect due to bed re-arrangement. However, detailed experiments by [START_REF] Shvidchenko | Flume study of the effect of relative depth on the incipient motion of coarse uniform sediments[END_REF] with well-sorted materials indicate that bed arrangement cannot entirely explain the increase in 𝜃 𝑐𝑟 . Indeed, the slope effect can be associated with changes in relative roughness 𝑘 𝑠 /ℎ, (with 𝑘 𝑠 the bed roughness height and ℎ the flow depth),

i.e., 𝑘 𝑠 /ℎ increases with 𝑆 [START_REF] Lamb | Is the critical shields stress for incipient sediment motion dependent on channel-bed slope?[END_REF][START_REF] Recking | Theoretical development on the effects of changing flow hydraulics on incipient bedload motion[END_REF][START_REF] Camenen | Discussion of "understanding the influence of slope on the threshold of coarse grain motion: Revisiting critical stream power[END_REF]. In a larger extent, the hydrograph can be related to the bed slope (steeper for high slopes) and flow acceleration may have an impact on friction, and thereby on bed shear stress [START_REF] Camenen | Discussion of "measurements of sheet flow transport in acceleration-skewed oscillatory flow and comparison with practical formulations[END_REF]. However, only steady flows will be considered in our study, in which case the slope effect can be regarded as a combination of at least the following two factors: bed arrangement and relative roughness.

As the direct parameters describing bed arrangement are often not reported in previous studies, we explore bed arrangement only through the effect of bed slope on the critical shear stress.

This study aims at discussing the estimation of the critical Shields number 𝜃 𝑐𝑟 and associated uncertainty. The paper is organised as follows: Section 2 is a review of existing methodologies for computing bed shear stress and critical value for inception of motion. Data collection and three 𝜃 𝑐𝑟 -predictive models based on 𝑑 * , 𝑆 or both are presented in Section 3 together with the Bayesian framework for uncertainty quantification. In Section 4, model parameters are estimated through Bayesian inference using prior knowledge on those parameters and observational data. The final estimation of 𝜃 𝑐𝑟 is then associated with a parametric uncertainty (related to the parameter estimation) and with a structural uncertainty (related to the choice of the model), which enables the evaluation of performance of the models. Results are discussed in Section 5, followed by 4 Perret, September 2, 2022 concluding remarks in Section 6.

SOURCES OF UNCERTAINTY IN MEASUREMENT OF CRITICAL BED SHEAR STRESS

Methods for bed shear stress computation

Various methods are available to compute the bed shear stress 𝜏 and most of them are reported in [START_REF] Wilcock | Estimating local bed shear stress from velocity observations[END_REF][START_REF] Biron | Comparing different methods of bed shear stress estimates in simple and complex flow fields[END_REF]. Different 𝜏 values (Shields number 𝜃 values) probably can be obtained depending on the chosen method. Those differences may explain a part of the scatter in the data of the critical Shields number.

Field studies demonstrate that the bed shear stress calculated from the depth-slope equation is generally larger than the one computed from the analysis of velocity profile [START_REF] Petit | Evaluation des critères de mise en mouvement et de transport de la charge de fond en milieu naturel[END_REF]). The first method provides a value at the cross-sectional scale 𝜏 𝑡 that lumps several components of flow friction such as the grain resistance 𝜏 ′ , which is responsible for inception of motion and bedload transport, and the bedform resistance 𝜏 ′′ (i.e., 𝜏 𝑡 = 𝜏 ′ + 𝜏 ′′ ). On the contrary, the velocity profile method yields the local bed shear stress, which can be assimilated to 𝜏 ′ . In most existing studies, indication about bedforms are almost missing; 𝜏 ′′ remains difficult to estimate and can represent 10 -75 % of 𝜏 𝑡 [START_REF] Buffington | A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers[END_REF]. According to [START_REF] Petit | Critical specific stream power in gravel-bed rivers[END_REF], an uncertainty of 50 % can be obtained for 𝜃 𝑐𝑟 if the calculation is based on the total bed shear stress

𝜏 𝑡 .
The major source of uncertainty for the depth-slope equation is mainly due to the estimation of the energy slope. For laboratory cases, the flume can even be too short to observe a water elevation gradient larger than the precision of the measuring device. In field cases, the depth-slope method is often improperly used, leading to large uncertainty, e.g., when the flow is not uniform or by replacing 𝐽 and 𝑅 ℎ by 𝑆 and ℎ, respectively. 5 Perret, September 2, 2022

One of the main difficulties using local methods in small scale laboratory experiments is to define the flow depth ℎ related to the reference bed level 𝑧 𝑏 , especially for coarse sediments for which a spatial variability does exist even if the bed is flat. [START_REF] Wilcock | Estimating local bed shear stress from velocity observations[END_REF] found that measurement uncertainty related to the velocity profile analysis and friction law methods was 5 % and between 5 -15 %, respectively. [START_REF] Biron | Comparing different methods of bed shear stress estimates in simple and complex flow fields[END_REF] ranked the Reynolds stress analysis as the most accurate method for beds with no forms and specific grain arrangement. For complex beds, the turbulent kinetic energy method (TKE) was recommended [START_REF] Kim | Estimating bottom stress in tidal boundary layer from acoustic Doppler velocimeter data[END_REF]. For velocity profile and friction law methods, one major issue is the definition of the roughness length 𝑍 0 .

Definition of incipient motion

One of the main issues related to bedload is the definition of incipient motion. Some exhaustive reviews [START_REF] Lavelle | Do critical stress for incipient motion and erosion really exist?[END_REF][START_REF] Dey | Sediment threshold[END_REF][START_REF] Beheshti | Analysis of threshold and incipient conditions for sediment movement[END_REF] for laboratory experiments, but applied their own definition for the incipient motion [START_REF] Kramer | Sand mixtures and sand movement in fluvial models[END_REF][START_REF] Vanoni | Measurements of critical shear stress for entraining fine sediment in a boundary layer. pasadena, california[END_REF]. Conducting experiments on a mixture of poorly sorted sand, [START_REF] Kramer | Sand mixtures and sand movement in fluvial models[END_REF] proposed the following four levels of sediment transport: (i) No transport, (ii) Weak transportfew of the smallest particles are in motion at isolated spots, (iii) Medium transport -particles of mean diameter are in motion at a small rate; and (iv) General transport -all particles are moving at all spots and at all times over the bed. [START_REF] Recking | Feedback between bed load transport and flow resistance in gravel and cobble bed rivers[END_REF] merged the second and third levels. [START_REF] Kramer | Sand mixtures and sand movement in fluvial models[END_REF] defined the threshold of motion to be the bed shear stress yielding general transport.

The main difficulty of the visual method is the distinction between the above levels. [START_REF] Vanoni | Measurements of critical shear stress for entraining fine sediment in a boundary layer. pasadena, california[END_REF] defined the threshold of incipient motion as the condition under which at least one grain is in movement every two seconds at any location. [START_REF] Neill | Quantitative definition of beginning of bed movement[END_REF] proposed a similar definition based on a dimensionless parameter 𝜖 = (𝑛Δ𝑡/𝐴) [𝜌𝑑 5 /(𝜌 𝑠 -𝜌)𝑔] 1/2 , where 𝑛 is the number of moving particles during a given time of observation Δ𝑡 on an observed bed area 𝐴. According to [START_REF] Neill | Quantitative definition of beginning of bed movement[END_REF], 𝜖 = 10 -6 corresponds to the inception of movement (≈ 0.8 grain/m 2 /s).

One issue remains: the validity of such criteria for any grain size.

These concepts of sediment threshold leads to a large scatter in the dataset and make comparisons difficult. It is obvious that there is no equivalence between the existing definitions. For example, both [START_REF] Vanoni | Measurements of critical shear stress for entraining fine sediment in a boundary layer. pasadena, california[END_REF]

(𝑞 𝑠 * ,ref ≈ 2.5 × 10 -4 ).

Evaluation of uncertainty in measurements of critical Shields number

It is possible to attribute an estimation of uncertainty to each data point 𝜃 𝑐𝑟 according to two uncertainty sources: definition of threshold for sediment motion (Δ def ) and methodology for computing bed shear stress (Δ 𝜏 ). The final uncertainty on 𝜃 𝑐𝑟 can be written as follows: As explained in Section 2, using a reference transport rate as an incipient motion definition is more robust than a visual definition (Δ def1 < Δ def2 , subscripts def1 and def2 correspond to reference transport rate definition and visual definition, respectively). Using the reference transport rate, the uncertainty lies mainly in the arbitrary chosen value for the reference transport rate but also in the reliability of measurements. For example, data collected with a Helley-Smith sampler and averaged throughout the river cross-section can lead to significant uncertainties [START_REF] Venditti | Bed load bias: Comparison of measurements obtained using two (76 and 152 mm) Helley-Smith samplers in a gravel bed river[END_REF][START_REF] Liu | Measuring bedload in gravel-bed mountain rivers: averaging> methods and sampling strategies[END_REF]. For laboratory experiments, bedload transport is often measured using a (2022), who observed that the visual definition generally leads to lower critical bed shear stress compared to other methods. Visual definition is generally not used for field data. There exists a data set from [START_REF] Young | Erosion velocities of skeletal carbonate sands, St. Thomas, Virgin Islands[END_REF] for which inception of motion was revealed by photo analysis.

Δ𝜃 𝑐𝑟 = 𝑢 𝜃 𝑐𝑟 𝜃 𝑐𝑟 = √ Δ 2 def + Δ 2 𝜏 (2)
We set the uncertainty for this case at 30%.

8

Perret, September 2, 2022

Uncertainty associated to the depth-slope method for bed shear stress computation is set for field experiments at 25 %, as the bed slope is generally used instead of the free surface slope and water depth may vary significantly throughout the river cross-section. These values are lower for laboratory studies, where bed conditions are constrained by the flume. The uncertainty for flume study Δ 𝜏,𝐷𝑆 is mainly linked to the calculation of the energy slope; we suggest Δ 𝜏,𝐷𝑆 =15%. This uncertainty may increase for specific cases with relatively coarse sediments and low water depths for which the spatial variability of the water depth is higher.

MATERIAL AND METHODS

Data compilation

We compiled an up-to-date data set for the estimation of 𝜃 𝑐𝑟 . It includes the data collected by [START_REF] Buffington | A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers[END_REF] but excluding data that used the competence function or theoretical developments. Indeed, the latter have not really been validated and lead to substantially different results compared to those obtained by the other methods. It also includes additional data collected by [START_REF] Recking | Theoretical development on the effects of changing flow hydraulics on incipient bedload motion[END_REF] as well as some additional data from the following studies : [START_REF] Rao | Stability and mobility of sand-bed channels affected by seepage[END_REF], [START_REF] Gregoretti | The initiation of debris flow at high slopes: experimental results[END_REF], [START_REF] Shvidchenko | Critical shear stress for incipient motion of sand/gravel streambeds[END_REF], [START_REF] Pilotti | Beginning of sediment transport of incoherent grains in shallowshear flows[END_REF] variability of a natural river that can also be biased due to scale effects, they remain of interest to study bedload transport processes since they provide data for controlled conditions with reduced uncertainties.

The measurement uncertainties Δ𝜃 𝑐𝑟 were estimated for each data point based on values reported in Table 2. Consequently, we obtained Δ𝜃 𝑐𝑟 = 21% for laboratory data using the reference transport rate definition, Δ𝜃 𝑐𝑟 = 25% for laboratory data using the the visual definition, and Δ𝜃 𝑐𝑟 = 35%

for field data (using the the reference transport rate definition). Since most of the data are from laboratory experiments, we eventually have an averaged value Δ𝜃 𝑐𝑟 = 25 %. Figure 1 presents the uncertainty 𝑢 𝜃 𝑐𝑟 as a function of the dimensionless grain size 𝑑 * . The largest values are observed for very fine and very coarse sediments since 𝜃 𝑐𝑟 can be over 0.1 for these specific grain sizes.

Tab. 2 here.

Fig. 1 here.

Again, this evaluation of the measurement uncertainties in data corresponds to a first rough estimation. The impact of the choice for Δ𝜃 𝑐𝑟 on the results is discussed in Section 5 using a sensitivity analysis.

Models for estimating 𝜃 𝑐𝑟

We propose here to test simple models for the estimation of the critical Shields number for inception of movement. First, we assumed the Shields curve can be evaluated as a function of the grain size only (through the input parameter 𝑑 * ) based on the equation of [START_REF] Soulsby | Threshold of sediment motion in coastal environment[END_REF]:

𝜃 (1) 𝑐𝑟 = 𝜈 (1) 1 𝑑 * + 𝜈 (1) 2 1 -exp(𝜈 (1) 3 𝑑 * ) (3)
where 𝜃 (𝑘) 𝑐𝑟 is the critical Shields number predicted by the model 𝑘 (here 𝑘 = 1), 𝜈 (1) 1 , 𝜈 (1) 2 , and 𝜈 (1) 3 are the parameters to evaluate (𝜈 (1) 1 = 0.24, 𝜈 (1) 2 = 0.055, and 𝜈 (1) 3 = -0.02 according to Soulsby 10 Perret, September 2, 2022

and Whitehouse (1997)). The [START_REF] Soulsby | Threshold of sediment motion in coastal environment[END_REF] equation was chosen since it is a continuous, single equation suitable for all grain size while including three fitting parameters only.

As compared to other formulas describing the empirical Shields curve as a function of the grain size [START_REF] Iwagaki | Fundamental study on critical tractive force[END_REF][START_REF] Van Rijn | Sediment transport, part I : bed load transport[END_REF], the Soulsby and Whitehouse (1997) equation yields very similar results. Some difference appears for the extrapolation for very fine sediments for which there is a lack of data for non-cohesive sediments. We assume a critical bed shear stress independent of the grain size as proposed by [START_REF] Soulsby | Threshold of sediment motion in coastal environment[END_REF].

Following the same idea, the critical Shields parameter can be evaluated as a function of the bed slope only, based on Recking (2009) equation:

𝜃 (2) 𝑐𝑟 = 𝜈 (2) 1 𝑆 + 𝜈 (2) 2 (4)
where 𝜈 (2) 1 and 𝜈 (2) 2 are the parameters to evaluate (𝜈 (2) 1 = 0.3 and 𝜈 (2) 2 = 0.04 according to Recking ( 2009)).

We propose to use a combination of Eqs. 3 and 4 to evaluate the critical Shields parameter as a function of both grain size and slope:

𝜃 (3) 𝑐𝑟 = 𝜈 (3) 1 𝑆 + 𝜈 (3) 2 × 𝜈 (3) 3 𝑑 * + 𝜈 (3) 4 1 -exp(𝜈 (3) 5 𝑑 * ) (5)
where

𝜈 (3) 1 , 𝜈 (3) 2 , 𝜈 (3) 3 , 𝜈 (3)
4 , and 𝜈 (3) 5 are the parameters to evaluate. Eq. 5 is an adjustment of Eq. 3 with an additional slope parameter. Eq. 5 is close to the following equation proposed by [START_REF] Camenen | Discussion of "understanding the influence of slope on the threshold of coarse grain motion: Revisiting critical stream power[END_REF]:

𝜃 𝑐𝑟 = 0.5 + 6𝑆 0.75 sin(𝜙 𝑠 -arctan 𝑆) sin(𝜙 𝑠 ) 0.24 𝑑 * + 0.055 [1 -exp(-0.02𝑑 * )] (6)
where 𝜙 𝑠 is the angle of repose of sediment. It should be noted that Eq. 5 does not include the possible instability due to steep slopes as Eq. 6 does. However, our data set is limited to bed slopes below 30% , above which the term sin(𝜙 𝑠arctan 𝑆)/sin(𝜙 𝑠 ) starts to be significant.

11

Perret, September 2, 2022

Bayesian estimation of predictive models

Overview and inference setup

Several sources of uncertainty affect the use of models in Eqs. (3-5). First, their parameters

𝜈 (𝑘)
𝑖 are unknown and will remain uncertain even after model calibration (parametric uncertainty).

Second, model calibration makes use of observed 𝜃 𝑐𝑟 that are uncertain as described in section 3 (observation uncertainty). Finally, the models are not perfect and are not expected to exactly replicate 𝜃 𝑐𝑟 (structural uncertainty).

Bayesian estimation provides a general and rigorous mechanism to estimate the unknown parameters of a model. It combines the information brought by uncertain calibration data with any pre-existing 'prior' information on the parameters. The method used in this paper is presented in details by Le Coz et al. ( 2014), [START_REF] Mansanarez | Bayesian analysis of stage-fall-discharge rating curves and their uncertainties[END_REF][START_REF] Perret | A rating curve model accounting for cyclic stage-discharge shifts due to seasonal aquatic vegetation[END_REF]. It was initially implemented for hydrometric rating curves but it can be applied to any kind of models.

Let 𝑶 = 𝑑 * ,𝑖 , 𝑆 𝑖 , 𝜃 𝑐𝑟,𝑖 , Δ 𝜃 𝑐𝑟 ,𝑖 𝑖=1,...,𝑛 denotes the 𝑛 observations in the dataset described in section 3. Each observation vector comprises values for grain size, slope, critical bed shear stress and its uncertainty (as described in section 2).

In addition, let 𝑀 denotes any of the models proposed in equations (Eqs. 3-5) to estimate a critical Shields number 𝜃 𝑐𝑟 from grain size 𝑑 * and/or slope 𝑆, with parameters 𝝂: 12 Perret, September 2, 2022

𝜃 𝑐𝑟 = 𝑀 (𝑑 * , 𝑆; 𝝂) (7)

Error model

The following error model is used to link observed and predicted values of 𝜃 𝑐𝑟 :

𝜃 𝑐𝑟,𝑖 = 𝑀 𝑑 * ,𝑖 , 𝑆 𝑖 ; 𝝂 𝜃 𝑐𝑟 ,𝑖 +𝛿 𝑖 + 𝜀 𝑖 (8)
This equation describes two distinct error sources. The error 𝛿 𝑖 is a measurement error and is assumed to be a realization from a Gaussian distribution with zero mean and known standard deviation Δ 𝜃 𝑐𝑟 ,𝑖 as described in section 2. The error 𝜀 𝑖 is a structural error due to the imperfection of the model 𝑀. It is also assumed to be a realization from a Gaussian distribution with zero mean.

However, its standard deviation 𝜎 is unknown, and therefore needs to be estimated along with parameters 𝝂. The reason behind this distinct treatment of observation and structural errors is that the former exists independently of any model, and its properties can therefore be estimated beforehand.

By contrast, the structural error is relative to the model of interest, and it is therefore difficult to know its properties before model estimation. Note that the normality of both measurement and structural errors is an assumption that can be evaluated through parameter estimation by examining residuals (i.e., observed minus predicted values). This assumption was found to be adequate for the data and models analysed in this work (not shown).

Prior distributions

For each unknown parameter, prior knowledge is encoded in a Gaussian distribution N (𝑚, 𝑠).

The mean value 𝑚 represents a 'prior guess' and the standard deviation 𝑠 represents the uncertainty around this prior guess. This standard deviation could potentially be very large when little is known about the parameter. In this paper, we choose to use as prior guess the values proposed by [START_REF] Soulsby | Threshold of sediment motion in coastal environment[END_REF] and [START_REF] Recking | Theoretical development on the effects of changing flow hydraulics on incipient bedload motion[END_REF] for Eqs. 3 and 4. For Eq. 5, prior guess values were evaluated assuming Eq. 5 corresponds to an adjustment of Eq. 3 using the additional slope parameter; so the prior guess values for parameters related to grain size were chosen equal to those of Eq. 3. Depending on the sensitivity on each of these parameters, a standard deviation was given between 30 and 50%. All prior guess values and related standard deviation are presented in Tab. 3 13 Perret, September 2, 2022

Tab. 3 here.

Outcome of Bayesian estimation

The raw outcome of Bayesian estimation is the posterior distribution of unknown parameters (𝝂, 𝜎). The probability density function (pdf) of this posterior distribution can be computed as shown in Appendix II. However, the posterior pdf is multi-dimensional and is therefore not easy to manipulate. Instead, it is more convenient to simulate many values from the posterior distribution, representing the posterior uncertainty in parameters. This simulation can be achieved by means of a Markov Chain Monte Carlo (MCMC) sampling algorithm. The particular sampler used in this paper is described in details in [START_REF] Renard | An application of Bayesian analysis and Markov chain Monte Carlo methods to the estimation of a regional trend in annual maxima[END_REF].

Once many values 𝝂 𝑗 , 𝜎 𝑗 𝑗=1,...,𝑁 𝑠𝑖𝑚 have been simulated by MCMC, the uncertainty in critical bed shear stress can be quantified by propagating these simulated values through the model: this corresponds to the Monte Carlo propagation method described in uncertainty analysis standards [START_REF] Jcgm | Evaluation of measurement data -Supplement 1 to the "Guide to the expression of uncertainty in measurement" -Propagation of distributions using a Monte Carlo method[END_REF]. In particular, applying the model equation (Eq. 7) 𝑁 𝑠𝑖𝑚 times yields 𝑁 𝑠𝑖𝑚 values of 𝜃 𝑐𝑟 that represent parametric uncertainty, i.e., the uncertainty due to the imperfect estimation of parameters 𝝂. The total uncertainty is obtained by adding to each of these 𝑁 𝑠𝑖𝑚 values a structural error 𝜀 randomly sampled from a Gaussian distribution with zero mean and standard deviation 𝜎 𝑗 . Note that measurement errors are not propagated at this stage, since the objective is to estimate the true 𝜃 𝑐𝑟 , rather than an observed, error-affected one. However, measurement errors still play an indirect role by affecting the posterior distribution and hence the uncertainty in estimated parameters.

RESULTS: ANALYSIS OF CRITICAL SHIELDS NUMBER UNCERTAINTIES

Evaluation of total uncertainty on critical Shields number using grain size only

Figure 3 plots the estimated 𝜃 (1) 𝑐𝑟 and related uncertainties using Eq. 3. Uncertainty bars of each data point are not plotted for the sake of readability. The best fit for Eq. 3 is obtained with As shown in Fig. 3, the total uncertainty originates mainly from the structural error and can be enclosed between -60% and +60%. The chosen model is certainly not the most appropriate one, i.e., 𝜃 𝑐𝑟 is not only function of the parameter 𝑑 * . A more appropriate model would yield a dominance of parametric errors, meaning that the uncertainty comes mainly from data.

𝜈 (1) 1 = 0.196, 𝜈 (

Evaluation of total uncertainty on critical shields number using bed slope only

Figure 4 presents the results obtained for the 𝜃 (2) 𝑐𝑟 = 𝑓 (𝑆) relationship and related uncertainties.

Here, the best fit for Eq. 4 is obtained with 𝜈 (2) 1 = 0.327 and 𝜈 (2) 2 = 0.0352, which is quite close to the results from Recking (2009) (𝜈 (2) 1 = 0.3 and 𝜈 (2) 2 = 0.04). Again, the total uncertainty comes mainly from the structural error and can be enclosed between -55% and +50%. For steep slopes (𝑆 > 0.1), the total uncertainty is lower and can be enclosed between -25% and +30%; the parametric error is no more negligible, meaning the model is more accurate here.

Fig. 4 here.

Evaluation of total uncertainty on critical shields number using both grain size and bed slope

Figure 5 depicts the results obtained for 𝜃 (3) 𝑐𝑟 = 𝑓 (𝑑, 𝑆) relationship and related uncertainties.

The best fit for Eq. 5 is obtained with

𝜈 (3) 1 = 1.055, 𝜈 (3) 2 = 0.274, 𝜈 (3) 3 = 0.510, 𝜈 (3
) 4 = 0.134, and 𝜈 (3) 5 = -0.068. These values are quite different to our prior guess, but this is not surprising, since we assumed the slope to be an adjustment coefficient of the critical bed shear stress evaluated as a function of 𝑑 * . Interestingly, the curve for the range 3 < 𝑑 * < 40 (i.e., sand-sized particles) is smoothed; the impact of grain size on the critical Shields parameter appears to be simpler than estimated from the Shields curve, i.e., inversely proportional to 𝑑 * for 𝑑 * < 3, and independent of 𝑑 * for 𝑑 * > 40. Indeed, most of data with sand particles were collected in low slope channels whereas those with gravel particles correspond to larger slopes. A fit without accounting for the slope effects is thus biased by the data collection.

Fig. 5 here.

In Fig. 5 are presented results for four specific slopes: 𝑆 = 0.001, 𝑆 = 0.02, 𝑆 = 0.1, and 𝑆 = 0.2. The plotted experimental data correspond to slope values of the same order (±25%); they are plotted with their uncertainties. Compared to Fig. 3, the total uncertainty is slightly reduced and can be enclosed between -50% and +50%. The total uncertainty is still dominated by the structural error. However, the parametric error becomes less negligible for highest slopes. It should be noted, however, that the proposed model underestimates 𝜃 𝑐𝑟 values for steep bed slopes, which could be due to the relatively low number of data describing high slopes.

Performance of models

Table 4 presents the statistical results of the predictive capabilities of the different equations (see graphs in Appendix III). 𝐸 𝑟,20 and 𝐸 𝑟,50 correspond to a percentage of data predicted accurately with allowed error of a factor 1.2 and 1.5, respectively; mean 𝑙𝑜𝑔 and std 𝑙𝑜𝑔 correspond to the mean and standard deviation of the logarithm of the ratio between the predicted and measured value.

It can be observed that the Bayesian inference leads to better predictive performances for Eqs. 3 and 4 since calibrated to the present larger data set as compared to the original calibration. A formulation with five parameters (Eq. 5) does not significantly improve the results apart for the standard deviation. When comparing measured to predicted 𝜃 𝑐𝑟 -values (see Figs. 8 and9), one can observe that all formulas (apart from the Camenen (2012) formula, which presents a larger scatter) yield relatively constant values, while observations vary a lot. This suggests that grain size and slope are not the only parameters to consider for predicting the inception of transport.

Tab. 4 here.

It is interesting to note that statistics presented in Tab. 4 slightly differ if we consider laboratory data or field data only (see also Tab. 5). The dispersion is higher for field data than for laboratory 16 Perret, September 2, 2022 data, as expected but the mean 𝑙𝑜𝑔 are also higher suggesting that 𝜃 𝑐𝑟 -values are found smaller in the field. Field data correspond generally to poorly sorted sediments. The median grain size may be not adapted or sufficient to characterize the inception of motion of the mixture. [START_REF] Recking | Theoretical development on the effects of changing flow hydraulics on incipient bedload motion[END_REF] suggested using 𝑑 84 instead of 𝑑 50 for poorly sorted sediments. However, the fine sediment fraction may be the key parameter for reducing the critical bed shear stress [START_REF] Wilcock | Methods for estimating the critical shear stress of individual fractions in mixed-size sediment[END_REF]. In a similar way, the visual observation definition yields in general smaller 𝜃 𝑐𝑟 -values than the reference transport rate definition [START_REF] Vah | A visual method for threshold detection of sediment motion in a flume experiment without human interference[END_REF]. This confirms the discussion in Section 2 and our suggestion to use a larger uncertainty for these data.

DISCUSSION

Assessment of retained data uncertainty

In Section 2, we attempted to evaluate uncertainty related to the technical sources for the critical Shields number data set. However, values in Table 2 remain partially subjective and arguable. A sensitivity analysis is performed to identify the impact of the choice of the uncertainty values on the results. We therefore explore the results obtained for a reduced (by a factor 1.3 and a factor 2)

or an increased (by a factor 1.3 and a factor 2) data uncertainty. These changes would lead to an average data uncertainty Δ𝜃 𝑐𝑟 of 12 %, 19 %, 32 % or 50 %, respectively. Let's remind that the averaged uncertainty initially was estimated at 25 %.

Fig. 6 shows the variation on the total and parametric uncertainties (𝐸 𝑡𝑜𝑡 and 𝐸 𝑝𝑎𝑟 , respectively) evaluated as an averaged of the ratio between the envelopes 97.5 % and 2.5 % for each 𝑑-values.

Consequently, an absence of uncertainty would yield the value of 1. The total uncertainty clearly decreases with an increase in the data uncertainties since the latter explains a larger part of the residual scatter. On the other hand, a minima for the parametric uncertainty is observed for data uncertainties between 19 and 25 %. With a lower data uncertainty, the models are too restricted and unable to properly fit the data. With a larger data uncertainty, the estimated parameters of the models are too uncertain. This uncertainty value between 19 and 25% corresponds to an optimum to evaluate the parameters of the model. This is consistent with our first evaluation of the data average uncertainty (i.e., Δ𝜃 𝑐𝑟 = 25 %).

Fig. 6 here.

For each equation (Eqs. 3, Eq. 4, and Eq. 5), we also tested the impact of the data average uncertainty Δ𝜃 𝑐𝑟 on the model coefficient estimation (see Fig. 10). Results are highly sensitive to the specified data uncertainty. In particular, for Eq. 5, our choices of priors significantly affect results for large data uncertainties. Indeed, we assumed the slope term (𝜈 (3) 1 𝑆 + 𝜈 (3) 2 ) as a correction of Eq. 3 (with 𝜈 (3) 2 = 1 as an initial prior) whereas the Bayesian approach indicates something intermediate : the posterior distribution of 𝜈 (3) 2 is three times narrower than its prior whereas the posterior distributions of 𝜈 (3) 3 and 𝜈 (3) 4 are three times wider than their prior. This would suggest that slope and grain size effects are competing each other.

A simple model for critical shields number using both grain size and slope

As the combination of slope and grain size effects reduce the smallest values observed for sand-sized particles, we propose to evaluate the critical Shields parameter as a function of both grain size and slope using the following four parameter equation:

𝜃 (4) 𝑐𝑟 = 𝜈 (4) 1 𝑆 + 𝜈 (4) 2 × 𝜈 (4) 3 𝑑 * + 𝜈 (4) 4 ( 9
)
where 𝜈 (4) 1 , 𝜈 (4) 2 , 𝜈 (4) 3 , and 𝜈 (4) 4 are parameters to evaluate. Since Eq. 5 is very similar to Eq. 9, we use the same priors, i.e., 𝜈 (4) 1 = 0.3, 𝜈 (4) 2 = 1, 𝜈 (4) 3 = 0.24, 𝜈 (4) 4 = 0.055. The best fit for Eq. 9 is obtained with 𝜈 (4) 1 = 1.158, 𝜈 (4) 2 = 0.180, 𝜈 (4) 3 = 0.410, 𝜈 (4) 4 = 0.195.

In Fig. 7 are presented results for four specific slopes: 𝑆 = 0.001, 𝑆 = 0.02, 𝑆 = 0.1, and 𝑆 = 0.2 in a similar way as in Fig. 5 for Eq. 5. Uncertainties are not improved compared to Fig. 5; they are enclosed between -55% and +55%. Nevertheless, when comparing to data, it clearly indicates there is no need of using such a complex empirical function of the grain size combined with the bed slope (Eq. 5). In addition, the simplified equation (Eq. 9) yields the best predictive performance compared to the other models with more than 70 % of the data predicted accurately with an allowed error of a factor 1.5 as shown in Table 4.

18 Perret, September 2, 2022 Fig. 7 here.

CONCLUSION

A series of equations was proposed to estimate the critical Shield number with an evaluation of its uncertainty. The models were derived from classical equations for the inception of movement relating the critical Shields number to grain size, or to longitudinal bed slope, or both. A Bayesian approach was used to estimate the model parameters using prior knowledge and observational data collected in literature. The Bayesian framework takes into account the measurement errors of the critical Shields numbers for the computation and gives two resulting uncertainty : a parametric (i.e., related to parameter estimation) and a structural (i.e., related to model itself) uncertainty.

The main sources of measurement errors were reported and discussed, especially those related to the definition of the inception of motion and to the method used to compute bed shear stress.

Measurement uncertainty was evaluated to 25 % in average for our data set. A sensitivity analysis was performed to discuss and verify this assumption by examining the impact of a reduced or increased measurement error on the results. Eventually, the proposed model (Eq. 9) improved results for estimating the critical bed shear stress for well-sorted sediments compared to existing models. However, for poorly-sorted sediments, one should use this model with the median grain size and apply additional laws for hiding and exposure effect.

A parametric uncertainty of approximately 10 % was found for 𝜃 𝑐𝑟 computed with models based on grain size or based on bed slope only. Total uncertainty was always larger than 50 %, which indicates significant structural uncertainty. A combination of both equations provided slightly better results. It also showed that smaller 𝜃 𝑐𝑟 values observed by [START_REF] Shields | Anwendung der Ahnlichkeits-mechanik und der turbulenzforschung auf die geshiebebewegung [application of similarity principles and turbulence research to bed-load movement[END_REF] for sand particles may be a bias linked to the combination of grain size and slope effects. Eventually, a model based on four parameters and assuming a continuous decrease in 𝜃 𝑐𝑟 with an increasing grain size yields the best results. Significant uncertainty remains; the parametric uncertainty being always smaller than the structural uncertainty. This indicates that the grain size and the bed slope are insufficient to describe the inception of movement. A more accurate estimation of the inception of motion should integrate 19 Perret, September 2, 2022 other factors, such as parameters describing bed composition (presence of fine sediments in a coarse matrix) and bed arrangement (including bed roughness, grain orientation and characteristic lengths of bed structures) as suggested by [START_REF] Perret | How does bed surface impact bedload transport over gravel-bed rivers?[END_REF]. However, such parameters were not measured in existing experiments, which limits the development of new predictive models.

In Tab. 1 are the main methods for computing bed shear stress presented. Figure 9 : Predicted 𝜃 𝑐𝑟,𝑒𝑞 versus measured 𝜃 𝑐𝑟,𝑒𝑥 𝑝 values of 𝜃 𝑐𝑟 using Eq. 6 (a), Eq. 5 (b), and Eq. 9 (c). Fig. 9. Predicted 𝜃 𝑐𝑟,𝑒𝑞 versus measured 𝜃 𝑐𝑟,𝑒𝑥 𝑝 values of 𝜃 𝑐𝑟 using Eq. 6 (a), Eq. 5 (b), and Eq. 9 (c).

  scale positioned at the downstream end of the flume[START_REF] Aguirre-Pe | Particle densimetric froude number for estimatingsediment transport[END_REF][START_REF] Perret | Transport of moderately sorted gravels at low bed shear stresses : impact of bed arrangement and fine sediment inltration[END_REF]). The uncertainty Δ def1 is evaluated thus equal to 10 %. For field experiments, since bedload transport is usually measured partially by sampling a finite number of points throughout the river section, we evaluate Δ def1 = 20%. For laboratory data, we propose an uncertainty using the visual definition Δ def2 = 15 % based on Perret (2017)'s experiments. This was recently confirmed byVah et al.

  ,[START_REF] Dey | Incipient motion of gravel and coal beds[END_REF],[START_REF] Dancey | Probability of individual grain movementand threshold condition[END_REF],[START_REF] Aguirre-Pe | Particle densimetric froude number for estimatingsediment transport[END_REF],[START_REF] Mueller | Variation in the reference shields stress for bed load transport in gravel-bed streams and rivers[END_REF], Hoffmans (2010),[START_REF] Prancevic | Unraveling bed slope from relative roughness in initial sediment motion[END_REF],[START_REF] Rousar | Incipient motion of coarse uniform gravel[END_REF],[START_REF] Perret | How does bed surface impact bedload transport over gravel-bed rivers?[END_REF]. The final data set includes 921 points (329 points obtained with the bedload extrapolation definition, and 592 points obtained wih the visual method). Most of these data are from laboratory experiments (867 points). The 54 points corresponding to field measurements were mostly obtained fromMueller et al's (2015) study using coupled measurements of flow and bed load transport in 45 gravel-bed streams and rivers in western North America. Even if flumes represent a small patch of the temporal and spatial 9Perret, September 2, 2022

  Figure 2.

Fig

  Fig. 2 here.

  Fig. 3 here.

Figure 1 :

 1 Figure 1 : Uncertainty 𝑢 𝜃 𝑐𝑟 as a function of the dimensionless grain size 𝑑 * (def1: reference transport rate, def2: visual definition).

Figure 2 :

 2 Figure 2: Diagram of the Bayesian model to evaluate the critical Shields number: inputs, parameters, observations, and outputs.

Figure 3 :

 3 Figure 3 : Critical Shields parameter 𝜃 𝑐𝑟 for inception motion estimated using Eq. 3 versus dimensionless grain size 𝑑 * ; Uncertainty envelopes for parametric (𝐸 𝑝𝑎𝑟 ) and total uncertainty (𝐸 𝑡𝑜𝑡 ) were defined as 95 % credibility intervals (𝐸 𝑡𝑜𝑡 =𝐸 𝑝𝑎𝑟 + 𝐸 𝑠𝑡𝑟𝑢𝑐 , with 𝐸 𝑠𝑡𝑟𝑢𝑐 the structural uncertainty).

Figure 4 :

 4 Figure 4 : Critical Shields parameter 𝜃 𝑐𝑟 for inception motion estimated using Eq. 4 versus longitudinal bed slope 𝑆; Uncertainty envelopes for parametric (𝐸 𝑝𝑎𝑟 ) and total uncertainty (𝐸 𝑡𝑜𝑡 )were defined as 95 % credibility intervals.

Figure 5 :

 5 Figure 5 : Critical Shields parameter 𝜃 𝑐𝑟 for inception motion estimated using Eq. 5 versus dimensionless grain size 𝑑 * for different values of slope: 𝑆 = 0.001 (a), 𝑆 = 0.02 (b), 𝑆 = 0.1 (c), and 𝑆 = 0.2 (d); Uncertainty envelopes for the parametric (𝐸 𝑝𝑎𝑟 ) and total uncertainty (𝐸 𝑡𝑜𝑡 ) were defined as 95 % credibility intervals.

Figure 6 :

 6 Figure 6 : Evaluation of the total and parametric uncertainties in the models as function of the averaged uncertainties Δ𝜃 𝑐𝑟 in experimental data (For θ(3) 𝑐𝑟 , (a): 𝑆 = 0.001, (b): 𝑆 = 0.02, (c): 𝑆 = 0.10, (d): 𝑆 = 0.20).

Figure 7 :

 7 Figure 7 : Critical Shields parameter 𝜃 𝑐𝑟 for inception motion estimated using Eq. 9 versus dimensionless grain size 𝑑 * for different values of slope: 𝑆 = 0.001 (a), 𝑆 = 0.02 (b), 𝑆 = 0.1 (c), and 𝑆 = 0.2 (d); Uncertainty envelopes for parametric (𝐸 𝑝𝑎𝑟 ) and total uncertainty (𝐸 𝑡𝑜𝑡 ) were defined as 95 % credibility intervals.

Figure 8 :

 8 Figure 8 : Predicted 𝜃 𝑐𝑟,𝑒𝑞 versus measured 𝜃 𝑐𝑟,𝑒𝑥 𝑝 values of 𝜃 𝑐𝑟 using Eq. 3 with Soulsby and Whitehouse (1997) coefficients (a) or with the coefficient estimated with the Bayesian approach (b) and Eq. 4 with Recking (2009) coefficients (c) or with the coefficient estimated with the Bayesian

Figure 10 :

 10 Figure 10 : Boxplots of parameter estimations for Eqs. 3 (a), 4 (b), and 5 (c) depending on average uncertainties Δ𝜃 𝑐𝑟 in experimental data.

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig.4. Critical Shields parameter 𝜃 𝑐𝑟 for inception motion estimated using Eq. 4 versus longitudinal bed slope 𝑆; Uncertainty envelopes for parametric (𝐸 𝑝𝑎𝑟 ) and total uncertainty (𝐸 𝑡𝑜𝑡 ) were defined as 95 % credibility intervals.

Fig. 7 .

 7 Fig.7. Critical Shields parameter 𝜃 𝑐𝑟 for inception motion estimated using Eq. 9 versus dimensionless grain size 𝑑 * for different values of slope: 𝑆 = 0.001 (a), 𝑆 = 0.02 (b), 𝑆 = 0.1 (c), and 𝑆 = 0.2 (d); Uncertainty envelopes for parametric (𝐸 𝑝𝑎𝑟 ) and total uncertainty (𝐸 𝑡𝑜𝑡 ) were defined as 95 % credibility intervals.

  Table. 1 (see Supplementary material). The depth-slope equation for uniform flow yields a reachaveraged value for 𝜏. The 1D Barré Saint-Venant equation (BSV) is preferred for non-uniform and unsteady flows. Friction laws calculate 𝜏 from the depth-averaged velocity (locally measured) or the cross-sectional-averaged velocity. Local bed shear stresses can also be estimated based on velocity profile measurements using either time-averaged values or fluctuations

Table 2

 2 recaps the proposed uncertainties according to the type of data based on expertise and literature. Note that the focus should not be on uncertainty values but rather on how they can be compared to each other. Most values in Table2were evaluated during marginal flume tests carried out for[START_REF] Perret | Transport of moderately sorted gravels at low bed shear stresses : impact of bed arrangement and fine sediment inltration[END_REF] study where critical Shield numbers were estimated with the different methods and definitions. The other values were assumed based on literature review (see Sections 2 and 2).Uncertainty for field data is expected to be larger than for laboratory data, partly because in-situ measurements are more difficult to achieve, grain size distributions are poorly sorted and often spatially distributed, and cross-sections are irregular with possible bedforms. Also, since field

measurements are often achieved over large periods, the studied river section may encounter bed changes. It should be noted that most data based on field experiments used here are from

[START_REF] Mueller | Variation in the reference shields stress for bed load transport in gravel-bed streams and rivers[END_REF]

, who used

[START_REF] Parker | Bed load and size distribution in paved gravel-bed streams[END_REF]

's criteria and depth-slope method.

  Uncertainties for the other local methods vary according to the topographic complexity of the studied zone (see Section 2). Local

techniques require many measurement points to evaluate a spatial and time-averaged bed shear stress. The presence of poorly sorted sediments makes also difficult the evaluation of the roughness length and bed level. Such local measurements are less common in the field. The proposed values in Table

.

2 are based on Perret (2017)'s experiments.

  Existing methods for the computation of bed shear stress. . . . . . . . . . . . . . . 29 2 Estimation of each component of the uncertainty depending on the type of data. . . 36 3 Prior specifications of the empirical parameters for the models (𝑘 = 1 to 4) discussed in this paper. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4 Statistics on the proposed equations to estimate the critical Shields number. . . . . 38 5 Statistics on the proposed equations to estimate the critical Shields number for selected parts of the data set. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

	List of Tables
	1
	35	Perret, September 2, 2022

TABLE 2 .

 2 Estimation of each component of the uncertainty depending on the type of data. : reference transport rate, def2: visual definition, 𝐷𝑆: depth-slope, 𝐹 𝐿: friction law, 𝑉 𝑃: velocity profile analysis, 𝑇 𝑃: turbulent profile analysis, 𝐵𝑆𝑉: 1D Saint Venant equation

	Component Laboratory data	Field data
	Δ def	Δ def1 = 10 %	Δ def1 = 20 %
		Δ def1 = 20 %	Δ def2 = 30 %
	Δ 𝜏	Δ 𝜏,𝐷𝑆 = 15 % Δ 𝜏,𝐷𝑆 = 30 %
		Δ 𝜏,𝐹 𝐿 = 12 %	-
		Δ 𝜏,𝑉 𝑃 = 10 %	-
		Δ 𝜏,𝑇 𝑃 = 8 %	-
		Δ 𝜏,𝐵𝑆𝑉 = 10 %	-

def136

Perret, September 2, 2022

TABLE 3 .

 3 Prior specifications of the empirical parameters for the models (𝑘 = 1 to 4) discussed in this paper.

	Equations	(𝑚, 𝑠)-values for each parameter
	𝜈 (𝑘) 1 𝑐𝑟 (Eq. 3) (0.24, 0.1) (0.055, 0.02) (-0.02, 0.01) 𝜈 (𝑘) 2 𝜈 (𝑘) 3 𝜃 (1)	𝜈 (𝑘) 4	𝜈 (𝑘) 5
	𝜃 (2) 𝑐𝑟 (Eq. 4) (0.3, 0.1)	(0.04, 0.02)
	𝜃 (3) 𝑐𝑟 (Eq. 5) (0.3, 0.1)	(1,0.5)	(0.24, 0.05) (0.055, 0.02) (-0.02, 0.01)
	𝜃 (4) 𝑐𝑟 (Eq. 9) (0.3, 0.1)	(1,0.5)	(0.24, 0.05) (0.055, 0.02)

TABLE 5 .

 5 Statistics on the proposed equations to estimate the critical Shields number for selected parts of the data set. Uncertainty 𝑢 𝜃 𝑐𝑟 as a function of the dimensionless grain size 𝑑 * (def1: reference transport rate, def2: visual definition). Critical Shields parameter 𝜃 𝑐𝑟 for inception motion estimated using Eq. 3 versus dimensionless grain size 𝑑 * ; Uncertainty envelopes for parametric (𝐸 𝑝𝑎𝑟 ) and total uncertainty (𝐸 𝑡𝑜𝑡 ) were defined as 95 % credibility intervals (𝐸 𝑡𝑜𝑡 =𝐸 𝑝𝑎𝑟 + 𝐸 𝑠𝑡𝑟𝑢𝑐 , with 𝐸 𝑠𝑡𝑟𝑢𝑐 the structural uncertainty).

	Equations	𝐸 𝑟,20 𝐸 𝑟,50 mean 𝑙𝑜𝑔 std 𝑙𝑜𝑔
	reference transport rate
	Eq. 3 (SW)	38.6 68.7	0.014	0.215
	Eq. 3 (Paper)	32.2 66.0 -0.069 0.211
	Eq. 4 (Rec)	47.1 75.7 -0.003 0.177
	Eq. 4 (Paper)	42.2 76.6 -0.045 0.176
	Eq. 6 (Cam)	26.1 56.5 -0.121 0.204
	Eq. 5 (Paper)	38.3 71.1 -0.069 0.186
	Eq. 9 (Paper)	45.3 74.2 -0.053 0.178
		visual observations
	Eq. 3 (SW)	23.6 54.5	0.059	0.247
	Eq. 3 (Paper)	33.6 63.7 -0.016 0.241
	Eq. 4 (Rec)	31.5 62.7 -0.014 0.268
	Eq. 4 (Paper)	36.4 64.9 -0.057 0.267
	Eq. 6 (Cam)	28.8 58.5 -0.087 0.235
	Eq. 5 (Paper)	37.4 67.3 -0.024 0.215
	Eq. 9 (Paper)	40.3 69.3 -0.003 0.199
		laboratory data
	Eq. 3 (SW)	29.1 60.5	0.035	0.233
	Eq. 3 (Paper)	33.8 65.1 -0.041 0.228
	Eq. 4 (Rec)	38.2 68.3 -0.017 0.236
	Eq. 4 (Paper)	39.5 69.7 -0.059 0.236
	Eq. 6 (Cam)	27.8 57.9 -0.104 0.224
	Eq. 5 (Paper)	38.4 69.5 -0.044 0.201
	Eq. 9 (Paper)	43.1 72.4 -0.025 0.186
		field data
	Eq. 3 (SW)	25.9 44.4	0.165	0.273
	Eq. 3 (Paper)	22.2 55.6	0.062	0.277
	Eq. 4 (Rec)	18.5 51.9	0.091	0.266
	Eq. 4 (Paper)	22.2 59.3	0.046	0.263
	Eq. 6 (Cam)	29.6 55.6 -0.026 0.226
	Eq. 5 (Paper)	27.8 55.6	0.032	0.270
	Eq. 9 (Paper)	25.9 50.0	0.044	0.280
	SW: Soulsby and Whitehouse (1997), Rec: Recking (2009), Cam: Camenen (2012), Paper: from
	this paper.
		39		Perret, September 2, 2022

Fig. 2. Diagram of the Bayesian model to evaluate the critical Shields number: inputs, parameters, observations, and outputs.

  Predicted 𝜃 𝑐𝑟,𝑒𝑞 versus measured 𝜃 𝑐𝑟,𝑒𝑥 𝑝 values of 𝜃 𝑐𝑟 using Eq. 3 with[START_REF] Soulsby | Threshold of sediment motion in coastal environment[END_REF] coefficients (a) or with the coefficient estimated with the Bayesian approach (b) and Eq. 4 with[START_REF] Recking | Theoretical development on the effects of changing flow hydraulics on incipient bedload motion[END_REF] coefficients (c) or with the coefficient estimated with the Bayesian approach (d).

		(a)			(b)
	10 0 10 0	Er 20 =29% Er 20 =27.9%		10 0 10 0	Er 20 =33.1% Er 20 =37.7%
		Er 50 =59.5% Er 50 =57.8%			Er 50 =64.5% Er 50 =68.7%
		mean log =0.04 mean log =-0.1			mean log =-0.03 mean log =-0.04
		std log =0.24 std log =0.22			std log =0.23 std log =0.21
	10 -1 10 -1			10 -1 10 -1
	10 -2 10 -2 10 -2 10 -2	10 -1 10 -1	10 0 10 0	10 -2 10 -2 10 -2 10 -2	10 -1 10 -1	10 0
		(c) (c)			(d)
	10 0 10 0	Er 20 =37.1% Er 20 =42.1%		10 0	Er 20 =38.5%
		Er 50 =67.4% Er 50 =71%			Er 50 =69.1%
		mean log =-0.01 mean log =-0.02			mean log =-0.05
		std log =0.24 std log =0.19			std log =0.24
	10 -1 10 -1			10 -1
	10 -2 10 -2 10 -2 10 -2	10 -1 10 -1	10 0 10 0	10 -2 10 -2	10 -1
	Fig. 8. 48		Perret, September 2, 2022

Perret, September 2, 2022

Acknowledgments

This study has been supported by INRAE (formerly Irstea), the Rhône-Alpes region through the CMIRA ExploraPro financial support (B. Camenen), and the French National Research Agency (ANR) under the grant ANR-18-CE01-0019-01 (DEAR project), and EDF (Ph.D thesis of E.

Perret). We would like to thank Alain Recking and Ladislav Rousar for providing their data sets.

Data availability statement

All data, models, or code that support the findings of this study are available from the corresponding author upon request.

APPENDIX II. EQUATIONS FOR BAYESIAN ESTIMATION

The posterior pdf of unknown parameters (𝝂, 𝜎) given the observed dataset 𝑶 can be computed prior pdf (19) where the symbol ∝ means 'is proportional to'. The likelihood results from the error model in equation ( 8) and can be computed as follows:

where 𝑓 N (𝑧; 𝑚, 𝑠) denotes the pdf of the normal distribution with mean 𝑚 and standard deviation 𝑠, evaluated at value 𝑧.

The prior pdf is computed as:

where 𝑓 U denotes the pdf of a uniform distribution between 𝑎 and 𝑏 and the 𝑓 N terms in the product denote the pdf of the normal distribution used as prior for each parameter, as described in section For Eq. 3 (Fig. 10a), the estimated model coefficients 𝜈 (1) 1 , 𝜈 (1) 2 , and 𝜈 (1) 3 decrease when data uncertainty increases, i.e., with a smaller sensitivity of Eq. 3 to grain size. For Eq. 4 (Fig. 10b), while the coefficient 𝜈 (2) 1 is independent of data uncertainty, 𝜈 (2) 2 decreases with larger uncertainty, indicating somehow a smaller impact of the less numerous data for high slopes. For Eq. 5 (Fig. 10c), one can observe a minimum and a maximum at our reference evaluation of the data uncertainty for

2 and 𝜈 (3) 4 , respectively. 27.9 57.8 -0.099 0.225 Eq. 5 (Paper)

) 5 = -0.068 37.7 68.7 -0.040 0.206

Eq. 9 (Paper, Section 5)