An extrapolation algorithm for estimating river bed grain size distributions across basins

Values representing grain size distributions of stream reaches are essential for estimating sediment transport at the reach scale. Various modeling frameworks exist that attempt to simulate reach-scale sediment transport across entire drainage basins to characterize sediment dynamics at a watershed scale. Such frameworks require estimates of grain size at each reach. Because obtaining direct measurements at this scale is impractical and logistically difficult, methods to estimate or extrapolate grain size measurements are needed, however, few currently exist. Here I present an extrapolation algorithm that uses one or more pebble counts to extrapolate full grain size distributions to each reach of a drainage network. In addition to the pebble count measurements, the tool requires a stream network geospatial feature class, attributed with values for reach-averaged slope and some consistent measure of relative flow magnitude (or a proxy for flow). I tested the tool in a set of sub-watersheds in the Bitterroot River basin of western Montana, US, with varying valley morphologies, and compared predictions to measurements at 16 sites. When using multiple measurements for calibration, mean absolute percent error averaged 5.8 % of the measured grain sizes in the phi scale. When using a single measurement for calibration, error averaged 8.4 %.

1 Introduction

A characterization of the bed surface grain size distribution (GSD) of rivers is essential for applying sediment transport formulae to estimate transport rates and total yield. In cases where the bed can be assumed to be equally mobile, a single characteristic size (e.g. median grain size, D₅₀) is often used to represent the bed material (e.g., Meyer-Peter and Müller, 1948). In contrast, to apply fractional transport equations (e.g., Parker, 1990; Wilcock and Crowe, 2003), a GSD is needed. GSD data can be tedious and time consuming to collect, however, and spatial coverage within a given area can be limited by accessibility to channels. Application is, therefore, generally limited in spatial extent unless GSDs can be estimated for locations without direct measurements.

Two general methods are typically used for quantifying GSDs: pebble counts and photo sieving. Pebble counts (often referred to as Wolman pebble counts after Wolman, 1954) are a simple method of characterizing bed GSDs by randomly gathering particles from the bed surface within a defined extent, and measuring the b-axis diameter of each. Sampling bias is common as larger particles are easier to see and are typically touched first as they protrude from the bed more, especially in poorly sorted mixtures (Leopold, 1970; Wohl et al., 1996). To avoid sampling bias, protocols have been recommended when performing pebble counts. For example, a “random walk” method is often used, with the sampler moving randomly across the defined area, identifying the particle at the point of the foot at each step and then retrieving it, rather than reaching down and grabbing the first particle touched. Kondolf (1997) recommends performing individual pebble counts for unique populations (geomorphic units) within a reach and taking the area-weighted average to characterize the whole reach. Another method uses a grid system, identifying the particle directly beneath each intersection and retrieving it as part of the count (Bunte and Abt, 2001). Typical pebble counts consist of 100 particles; fewer may be needed to characterize relatively homogenous beds, and more may be needed in poorly sorted beds with a wide GSD. Additionally, multiple pebble counts can be performed for a given location by a single person or multiple individuals, and the results used to quantify uncertainty in the values used (e.g., D₅₀). Despite the drawbacks of this method related to sampling bias, it is easy to implement, provides good characterizations of river beds, and is the most widely used method.

Photo sieving is a methodology in which image processing algorithms are used to identify individual grains in photographs and calculate their diameter (Ibbeken and Schleyer, 1986). The advantages of this method are that rather than a random sample, all identifiable particles can be used to generate the GSD, and it can be applied over large areas wherever photographs are available. These algorithms have been improved upon with advances in computer vision technology (e.g., Buscombe, 2008), including even broader application with the advent of drone-based surveys, which can produce orthoimagery of very large areas of river bed (e.g., Purinton and Bookhagen, 2019). Drawbacks of this method are that it can only be effectively applied where grains are exposed and easily identifiable in images, and that partially buried or obscured grains that are still part of the surface layer cannot be accurately measured. This methodology, therefore, cannot be applied to inundated or densely vegetated areas. Where applicable, photo sieving produces results comparable to pebble counts (Strom et al., 2010). Where it is appropriate to use, photo sieving can characterize grain size over larger spatial scales than pebble counts, however, like pebble counts, it is an impractical method for quantifying GSDs across entire drainage networks.

Several additional methods exist for characterizing bed GSDs. For very fine sediments, laser diffraction can be used to measure particle size (Blott et al., 2004; Eshel et al., 2004). For fine beds (silt to gravel), samples can be taken in the field that are separated into different size fractions in the lab using sieves or settling tubes (Komar and Cui, 1984). When available, very high-resolution topography of the bed (e.g., from terrestrial or mobile laser scanners or structure-from-motion photogrammetry) to estimate median grain sizes based on relationships between grain size and roughness (Pearson et al., 2017; Williams et al., 2020).

Controls on grain size include various factors relating to sediment production processes, as well as hydraulic sorting of bed material. For example, lithology exerts control on sediment delivery rates to channels, as well as the size of that sediment (Mueller et al., 2016). The size of individual grains produced from a rock layer is also a function of fracture spacing within the rock (Neely and DiBiase, 2020; Verdian et al., 2021). When sediment is delivered to channels from hillslope erosion processes, hydraulic sorting then exerts control on what sizes are present in the bed at a given reach (Snelder et al., 2011). As a result, variables often associated with sediment transport (i.e., channel gradient and flow magnitude) also correlate with bed GSD (Cohen et al., 2022).

Various sediment models seek to simulate transport across drainage networks in a spatially explicit manner (e.g., Czuba, 2018; Gilbert and Wilcox, 2020; Schmitt et al., 2016), and thus require a specification of grain size at each spatially distinct location in the modeling domain. These models have traditionally used relatively simple approaches, including applying the same distribution to the entire network, or correlating a median grain size with topographic and hydrologic information. A systematic method for estimating full distributions across drainage networks would thus be a valuable contribution for improving existing sediment transport and connectivity models.

Several approaches have been used to estimate grain size across large spatial extents. Jha et al. (2022) assessed D₅₀ at over 2400 locations from five different sources and used inverse distance weighted interpolation to estimate D₅₀ regionally between sample points. This provides broad regional patterns, but does not provide information in specific catchments or reaches. Abeshu et al. (2022) attributed D₅₀ information from 2577 United States Geologic Survey (USGS) gages to corresponding reaches in the National Hydrography Dataset (NHD) and used machine learning to develop a model based on the attributes available in the in NHDPlus dataset, finding that slope and distance between the reach and the stream's mouth were the most important channel characteristic predictors of grain size. Similarly, Snelder et al. (2011) used random forest models to assess data from 677 stream reaches in France with grain size measurements and found that the model incorporating variables related to hydraulic sorting (e.g., slope and precipitation) best predicted D₅₀.

Other approaches for estimating grain size across large spatial scales have used remote sensing-based methods. For example, Snyder et al. (2013) used digital elevation models to predict grain size. Their approach requires spatially distributed estimates of Manning's n roughness values, and is based upon the assumption that critical Shields stress (${\mathit{\tau }}_{\mathrm{c}}^{\ast }$) is 0.04 and that it is approximately equal to the Shields stress of a 2-year recurrence-interval flood. Similarly, models based on Shields stress (Buffington et al., 2004) and stream power (Gorman et al., 2011), calibrated empirically with both GIS and field data, have been used to develop predictive grain size models. However, all of the approaches described here estimate only the D₅₀, and no method currently exists that predicts full GSDs across drainage networks, despite the fact that such information would be useful for a variety of purposes.

In this paper, I introduce a method for extrapolating grain size distributions for coarse-bedded rivers (sand and coarser) throughout a drainage network based on one or more measurements (e.g., obtained using pebble counts or photo sieving). This method, similar to previous approaches, has the advantage of being calibrated based on actual measurements for a basin, while utilizing remote sensing data to obtain topographic data. It also predicts a full grain size distribution for each reach, rather than just D₅₀, allowing for application of fractional sediment transport functions. It also essentially estimates grain sizes as a function of hydraulics. Therefore, because hillslope sediment delivery is an additional control on grain sizes present within a reach (Hassan et al., 2019), this approach is best suited to alluvial river systems, systems with minimal coarse sediment inputs from hillslopes, or systems where the magnitude of hillslope sediment production is not highly spatially variable.

2 Methods

The grain size tool presented here can use a single measurement, or multiple measurements of GSD to extrapolate estimates of D₁₆ (the particle size for which 16 % of the bed is finer), D₅₀, and D₈₄ (the particle size for which 84 % of the bed is finer), along with an uncertainty range for each to every segment or reach of a geospatial drainage network layer. The tool then generates an estimate of full GSDs for each reach (i.e., the fraction of each half-phi size class in the bed; this information is stored in a “.json” file). Because relative size differences are more meaningful than absolute size differences in GSDs (i.e., a difference of 1 mm from a 1 mm to a 2 mm particle has much greater implications for transport than from a 150 to 151 mm particle), the logarithmic transformation of particle sizes from standard units (mm) to phi units is often used (Tanner, 1969). As such, I record GSDs in the phi scale and report errors in the phi scale in the results. The tool requires an input drainage network GIS feature class, for which GSDs are estimated and attributed to each feature. The layer should include attributes representing reach-averaged slope and some measure of flow or a proxy for flow magnitude (e.g., a two-year recurrence interval flow or drainage area; I used cumulative annual precipitation within the basin upstream of each reach), with a value for each feature.

2.1 Algorithm

Two different processes can be used, depending on whether one or multiple pebble count measurements are being used for calibration. For best results, multiple measurements spanning a range of stream sizes and gradients should be used. However, data availability is often limited, sometimes to a single measurement for a basin, so the tool also uses a novel algorithm to estimate GSDs in such circumstances. The algorithm steps are listed below, followed by more in-depth descriptions.

Both (single measurement and multiple measurements):

• 1.

  Calculate D₁₆, D₅₀, and D₈₄ and their uncertainty for measurement reach(es) from GSD measurements (pebble counts).

Single Measurement:

• 2.

  Find the flow depth required to transport the D₁₆, D₅₀, and D₈₄ at the measurement reach (critical depth).

• 3.

  Find hydraulic geometry coefficients associated with the critical depth.

• 4.

  Extrapolate critical depth to each drainage network segment based on a parameterized hydraulic geometry relationship.

• 5.

  Estimate D₁₆, D₅₀, and D₈₄ at each segment using the Shields equation and a chosen method for estimating critical Shields stress.

Multiple:

• 6.

  Calibrate functions for D₁₆, D₅₀, and D₈₄ that relate grain sizes to the product of slope and flow using all available measurements.

• 7.

  Estimate D₁₆, D₅₀, and D₈₄ at each segment using calibrated empirical function.

Both:

• 8.

  Fit skew-normal distribution parameters (location, scale, and shape) for each pebble count measurement. Take the average value of each parameter to use as a starting point in estimating the parameters for non-measurement reaches.

• 9.

  For each drainage network segment, fit new skew-normal distribution parameters using the parameters from the measurement reaches as a starting point. These parameters are adjusted until the error between D₁₆, D₅₀, and D₈₄ found from generating a new distribution is minimized compared against the actual values attributed to the segment in step 5 or 7.

• 10.

  Generate a distribution using the solved parameters, removing fractions finer than sand for the segment. Calculate the fraction of the resulting distribution within each half-phi size interval and store the data in the output .json file.

2.1.1 Calculating D₁₆, D₅₀, and D₈₄ from measurements and quantifying uncertainty

For each pebble count, the data is bootstrapped with replacement 100 times and sampled until the number of values is 80 % of the total number of measurement values. For each bootstrap, the 16th, 50th, and 84th percentile of the data is recorded, generating a distribution of each consisting of 100 values. The values for D₁₆, D₅₀, and D₈₄ are then set as the mean of these distributions with the uncertainty being characterized as plus or minus two standard deviations. Additionally, the reach slope and flow attribute are recorded (from the drainage network attributes) in order to parameterize relationships between the product of slope and flow and the different grain sizes.

2.1.2 Critical Depth

In the case where a pebble count at a single location is used, after the D₁₆, D₅₀, and D₈₄ are calculated, the depth required to mobilize each of these sizes (critical depth) is estimated using the Gilbert and Wilcox (2024) method for calculating critical Shields stress, along with the Shields function. Critical Shields stress for a given grain size D_i is calculated as

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\mathrm{c},i}^{\ast }=a{\left({\displaystyle \frac{D_i}{D_{\mathrm{50}}}}\right)}^{-\mathrm{0.67}}\end{array}$$

where the subscript i denotes a given grain size fraction, and the coefficient a varies as a function of the D₈₄ to D₅₀ ratio,

$$\begin{array}{}\text{(2a)}& {\displaystyle }a=\mathrm{0.029}\mathrm{if}{\displaystyle \frac{D_{\mathrm{84}}}{D_{\mathrm{50}}}}\le \mathrm{2}\text{(2b)}& {\displaystyle }a=\mathrm{0.043}\ln \left({\displaystyle \frac{D_{\mathrm{84}}}{D_{\mathrm{50}}}}\right)-\mathrm{0.0005}\mathrm{if}{\displaystyle \frac{D_{\mathrm{84}}}{D_{\mathrm{50}}}}>\mathrm{2}\end{array}$$

The critical Shields number is given by the function

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{\mathrm{c},i}^{\ast }={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{c},i}}{\left({\mathit{\rho }}_{\mathrm{s}}-\mathit{\rho }\right)gD}}\end{array}$$

where τ_c,i is the critical boundary shear stress at which the ith size fraction of the bed becomes mobile. Shear stress is typically estimated using the depth-slope product (which assumes steady, uniform flow),

$$\begin{array}{}\text{(4)}& \mathit{\tau }=\mathit{\rho }ghS.\end{array}$$

In equations 3 and 4, g is acceleration due to gravity (9.81 m s⁻²), ρ_s is the density of sediment (2650 kg m⁻³), ρ the density of water (1000 kg m⁻³), h is average flow depth (m), and S is reach averaged slope (m m⁻¹). τcan be solved for by multiplying the denominator of the Shields function by τ^∗, and critical depth for the size fraction i can then be estimated by rearranging the depth-slope product,

$$\begin{array}{}\text{(5)}& h_{\mathrm{c},i}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{c},i}}{\mathit{\rho }gS}}.\end{array}$$

2.1.3 Hydraulic Geometry Coefficients

In order to extrapolate critical depth values to each drainage network segment, hydraulic geometry coefficients are parameterized. Hydraulic geometry equations scale channel geometry with discharge (Leopold and Maddock, 1953). In downstream hydraulic geometry relationships, an exponent of 0.4 is often used for depth (Morel et al., 2020):

$$\begin{array}{}\text{(6)}& h=\mathit{\alpha }{Q}^{\mathrm{0.4}}\end{array}$$

where Q is discharge. In this case, Q can be substituted for drainage area or whatever other flow proxy attribute is being used. α is simply solved for in the measurement reach, and used for every other drainage network segment.

2.1.4 Estimating Grain Size (Single Measurement)

After the coefficient in Eq. (6) is calculated, critical depth for D₁₆, D₅₀, and D₈₄ is estimated for each reach, based on the reach's flow proxy attribute value. Shear stress at that depth is calculated using Eq. (4) with the calculated critical depth, and then the Shields formula (Eq. 3) is rearranged to solve for D.

$$\begin{array}{}\text{(7)}& D_i={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{c},i}}{\left({\mathit{\rho }}_{\mathrm{s}}-\mathit{\rho }\right)g{\mathit{\tau }}_{\mathrm{c},i}^{\ast }}}\end{array}$$

where the subscript i represents the specific grain size percentile (i.e., 16, 50, 84). Section 2.1.1 through 2.1.4 are summarized graphically in Fig. 1.

Figure 1 A graphical representation of the algorithm for predicting grain sizes at all reaches of a drainage network using measurements from a single reach. D₁₆, D₅₀, and D₈₄ are calculated from a pebble count (1) at a measurement reach, the critical depth for that grain diameter estimated (2), and hydraulic geometry coefficients for that depth calculated (3). That coefficient is used to extrapolate critical depth to the remainder of the drainage network (4), and used in a rearrange Shield's Equation to estimate the corresponding grain size at every reach of the drainage network (5).

2.1.5 Estimating Grain Size (Multiple Measurements)

When multiple pebble counts are used with the tool, three functions are parameterized to calculate D₁₆, D₅₀, and D₈₄ respectively. For each of the three grain size fractions, a least squares regression is used to fit a relationship predicting grain size as a function of the product of slope and flow (Fig. 2). The interaction of these two variables is largely able to explain variation in grain size. When reaches are steeper, they tend to have coarser beds, as flow competence is higher and larger particles can be mobilized. However, a steep headwater stream with little flow will have a much finer GSD than a higher-order stream with the same slope, where the combination of slope and flow result in much greater flow competence. The values of slope and flow for every other reach in the drainage network are then used to predict D₁₆, D₅₀, and D₈₄ using the fitted relationships.

Figure 2 Parameterized relationships (lines) for calculating D₁₆, D₅₀, and D₈₄ as a function of the product of slope and flow using four pebble counts at different locations within a single basin in the study area. In this application, I used precipitation within the upstream basin as the flow proxy, therefore the x axis is the product of these precipitation values and the reach-averaged slope.

In the case of both single and multiple measurements, as D₁₆, D₅₀, and D₈₄ are calculated for each reach, the value is stored in the attribute table along with a minimum and maximum value (based on the uncertainty previously calculated following Sect. 2.1.1). These values are subsequently used to calculate the full GSDs for the reach.

2.1.6 Estimating Full Grain Size Distributions

For each grain size measurement, the data is converted to phi values. Typically, phi value distributions for coarse-bedded channels are negatively skewed. The parameters of a best-fit skew-normal distribution (location, scale, and shape) are therefore calculated (Fig. 3a). In the case of multiple measurements, the parameters are found for each measurement, and the mean value of each parameter is used as a starting point for estimating parameters at each individual reach.

Figure 3 Pebble count measurements are used to fit a skew-normal distribution and calculate its parameters (a). Those parameters are then used to generate new distributions in non-measurement (prediction) reaches (b). The parameters are then adjusted until the errors between D₁₆, D₅₀, and D₈₄ from newly generated distributions and previously predicted values for D₁₆, D₅₀, and D₈₄ are minimized.

For each network reach, a distribution is generated using the parameters calculated from the measurement data. An optimization process is then implemented, varying the parameters until the errors between D₁₆, D₅₀, and D₈₄ of a distribution generated from the new parameters and the values of D₁₆, D₅₀, and D₈₄ calculated for that segment previously (Sect. 2.1.4 or 2.1.5) are minimized (Fig. 3b). These optimized parameters are used to generate a final distribution. Values finer than sand are removed (in this application to coarse bedded rivers). From the remaining values, the fraction within each half-phi size interval is calculated, and the values stored in a .json file.

3 Application

To test its effectiveness, I applied the algorithm to two collections of mountain river sub-basins within the larger Bitterroot River basin of Western Montana (Fig. 4). Within each basin, I performed pebble counts at several locations spatially dispersed across the watershed, capturing varying stream sizes and gradients. One basin in each collection was used to calibrate the model, after which it was applied to the remaining basins. The first collection consisted of basins with dendritic networks flowing through v-shaped valleys, where steep hillslopes physically connected to channels are prevalent throughout the watersheds (Woods, Rye, Sleeping Child, and Burnt Fork Creeks). The second collection was made up of streams within formerly glaciated mountain basins (Blodgett, Roaring Lion, and Lost Horse Creeks), with u-shaped valleys and wider valley bottoms, and greater variability in gradients of hillslopes adjacent to channels as a result of complex erosion patterns and valley geometries common in post-glacial landscapes (Hassan et al., 2019). Because of the greater uniformity in hillslope gradients adjacent to channels in first group, I infer generally high hillslope-channel connectivity, and refer to this group hereafter as the “uniform group.” In the second group, because hillslope gradients adjacent to channels were lower and more variable, I infer lower and more variable hillslope-channel connectivity, and refer to this group hereafter as the “variable group.” After running the model using multiple calibration measurements, I ran it for each basin using a single calibration measurement within the basin, allowing for comparison of results using the two methods.

Figure 4 Locations of grain size sampling within the Bitterroot River basin, and the location of the Bitterroot basin within Montana. Coordinates are in WGS84 UTM 11N. Imagery: © Google Maps, 2025 Airbus.

4 Results

For model results using multiple pebble count measurements for calibration (Table 1), the average discrepancy between predicted and measured size was 6 mm for D₁₆, 14 mm for D₅₀, and 32 mm for D₈₄, or mean absolute percent error (MEAP) of 26.1 %, 25.8 %, and 26.1 % respectively. In comparing predictions to measurements in the phi scale, this equates to MEAP of 7.6 % for D₁₆, 5.3 % for D₅₀, and 4.4 % for D₈₄. All combined, the MEAP for phi values was 5.8 % (26 % for mm). With the exception of D₈₄, which had similar error in both cases, errors were lower for the uniform group (6 %, 4.1 %, and 5 % MEAP for phi values) than for the variable group (9.1 %, 6.4 %, and 3.9 % MEAP for phi values).

Table 1 Grain size estimate results using multiple calibration measurements.

Within the uniform group, the geology of three of the four watersheds is primarily metamorphic (gneisses and schists). By contrast, the fourth basin is underlain by sedimentary and metasedimentary (e.g., quartzite) rocks. If the two sites with different geology than the others (BF and GC) are excluded, error is further reduced (5.1 %, 1.9 %, and 2.5 % MEAP for phi values). In the variable group, errors were the highest where GSD was overpredicted in a very low-gradient, unconfined reach.

Using a single calibration measurement (Table 2), errors were slightly higher than they were using multiple calibration measurements (11, 21, 48 mm; 13.1 %, 6.6 %, 5.7 %). The average error of all measurements was 8.4 %.

Table 2 Grain size estimate results using single calibration measurements.

One useful feature of the model is that, as previously mentioned it attributes geospatial data with the grain size predictions for D₁₆, D₅₀, and D₈₄. Figure 5 highlights how this allows for identification of spatial patterns in grain size throughout a basin. As qualitative validation of these spatial patterns, obvious differences in grain size distributions of different reaches based on photographs are captured in the model predictions (Fig. 6).

Figure 5 The drainage network for the Sleeping Child Creek basin, with each reach attributed by its predicted D₅₀ value. The network is overlain on a hillshade to show the topography hypothesized to result in uniformly high hill-slope channel connectivity.

Figure 6 Photographs along different reaches of Roaring Lion Creek with the predicted D₁₆, D₅₀, and D₈₄ values.

The .json format is simply a form of text file that is useful for storing information in a human readable format that is also interoperable with various programming languages. The output .json from the tool stores all of the grain size information produced by the tool for each reach. Specifically, this includes the calculated D₁₆, D₅₀, and D₈₄, as well as the fraction of each half-phi size interval in the full GSD.

5 Discussion and Conclusion

One of the limitations of the tool is that by using variables related to flow competence to predict grain sizes, it assumes that grain size distributions are hydrologically driven. In mountain drainage basins, as previously mentioned, grain sizes can be driven by, and vary greatly as a function of hillslope-channel connectivity and associated erosion and transport processes. As a result, mountain drainage networks often include channels that are completely non-alluvial, whose bed composition is entirely a function of non-fluvial processes (e.g., glacial erosion, rockfall, or other hillslope sediment production processes). In these cases, where the GSD is entirely a function of hillslope sediment delivery with no fluvial reworking, this model is not applicable for predicting grain sizes. In our study basins, the majority of the sediment delivered from hillslopes to channels in the uniform group watersheds is mobile at some flood flows. As a result, our prediction errors were fairly uniform across all sites. In contrast, much of the sediment delivered from hillslopes to channels in the variable group basins is large enough to be immobile at all flows. As a result, prediction errors in these basins were more variable.

Overall, results from the application of this algorithm presented here suggest that even in mountainous watersheds where GSDs are also a function of hillslope processes, using a hydraulics-based approach for estimating GSDs provides reasonable estimates.

Estimates of GSDs produced using this method may vary based on the input data used and the method of data collection. Grain size values are a function of flow or a flow proxy. Because these empirical relationships are calculated as a function of these data, the specific metric used should not have a large effect on the outputs as long as the data is spatially consistent with similar accuracy across basins. The method used to collect GSD measurements may result in a different distribution shape, which would affect the distribution parameters used to generate GSDs in the prediction reaches. However, results using different methods are typically fairly similar (Strom et al., 2010). In the single measurement method, because a single hydraulic geometry coefficient is calculated from the measurement reach, if the reach is notably different from other reaches in the basin (e.g., much rougher, which could result in greater depth for the discharge value), the coefficient may be a poor fit for much of the basin and result in systematic under or overestimation of grain sizes in the rest of the reaches.

The highest prediction error for these basins occurred in a very low gradient reach where grain sizes were overestimated. This reach, although it contained hillslope derived boulders along channel margins, had a sand and pea gravel bed. This approach is designed for gravel bedded reaches where GSDs typically follow a skew-normal distribution, and this reach was an outlier within the watershed that likely fell outside of these assumptions. The approach could be modified (e.g., use a different distribution for different types of beds), however this highlights the importance of accounting for the assumptions of the tool when interpreting the outputs. In our results, there was no heteroscedasticity in the errors for D₁₆, D₅₀, and D₈₄ calculated as phi values (p = 0.13, Breusch-Pagan test), indicating that uncertainty is similar across the full GSD using this approach.

In this application, I used single pebble counts to represent the GSD of entire reaches in the interest of gathering broad spatial coverage. Reaches are better characterized by performing pebble counts within each distinct geomorphic unit present, and subsequently taking the area-weighted mean (Kondolf, 1997). In applications where higher accuracy is desired for a smaller geographic area, calibrating the model with these types of pebble counts may produce better results.

Comparing measured GSDs to aerial imagery for very coarse-bedded reaches suggests that I may have had a bias for the wetted channel at the time of sampling. Consequently, I may have underestimated the size distribution by oversampling the more mobile fraction present near the thalweg relative to other, coarser grains distributed throughout the channel (Fig. 7). In very coarse beds, where it is more challenging to walk over the whole bed, a random walk sampling frame is likely less applicable and less effective. In these cases, a different sampling frame (e.g. gridded sampling) could improve sampling of all sizes present in the bed.

Figure 7 Aerial imagery of a sample reach of Blodgett Creek with a scale bar. The wetted portions of the channel have a finer GSD than the bed as a whole.

Within the uniform group, when using a basin with metamorphic lithology for calibration, errors were the greatest in the one sub-basin underlain by different geology. This highlights the control geology plays on GSDs, and consequently sediment yield, in line with previous findings that geology is a primary control on sediment supply to channels in the Northern Rocky Mountains (Mueller and Pitlick, 2013). While geology exerts control at the landscape scale, other factors can influence variation and GSD at local scales; for example confinement (O'Brien et al., 2019), large wood (Welling et al., 2021; Wohl and Scott, 2017), and beaver dam building activity (Levine and Meyer, 2014). These factors are not included in the model, but should be considered when applying it.

The spatial variability in channel bed grain sizes arises from the interplay between sediment delivery, both from upstream and from adjacent hillslopes, and the competence of the flow to transport that material. These differences can create large variability in transport rates from reach to reach, driving heterogeneity in fluvial sediment connectivity with implications for landscape evolution, habitat, and human infrastructure. The algorithm presented here provides an easily applicable method to use pebble counts to extrapolate estimates for full GSDs to an entire drainage network, allowing for increased accuracy in watershed-scale sediment transport estimates.