Grain Size Distribution Extrapolation Algorithm

Introduction

The grain size tool uses a single or multiple measurements of grain size distribution (e.g., Wolman (1954) pebble count) to extrapolate estimates of $D_{16}$, $D_{50}$, and $D_{84}$ along with an uncertainty range for each to every segment or reach of a drainage network GIS layer. Additionally, the tool produces a .json file where full grain size distribution information is stored for each reach (i.e., the fraction of each size class in the bed). The steps are outlined below, followed by a more in-depth description of each.

Algorithm Steps

1. Calculate $D_{16}$, $D_{50}$, and $D_{84}$ and their uncertainty from grain size distribution measurements.
2. Find $D_{max}$ for measurements and parameterize a linear equation or scalar for estimating $D_{max}$.
3. Find the flow depth required to transport the $D_{16}$, $D_{50}$, and $D_{84}$ at measurement reaches.
4. Find hydraulic geometry coefficients associated with the critical depths.
5. Extrapolate critical depths to each drainage network segment based on parameterized hydraulic geometry relationship.
6. Estimate $D_{16}$, $D_{50}, $D_{84}$ at each segment using the Shields equation and the Gilbert (2022) method for estimate critical Shields stress. Estimate $D_{max}$ using the relationship developed in step 2.
7. Take the square root of the measurement data to normalize the original distributions and fit normal distribution parameters for each measurement. Take the average of each for the final parameters.
8. For each drainage network segment, fit new normal distribution parameters using the parameters from the measurement reaches as a starting point. These parameters are adjusted until the error between $D_{16}$, $D_{50}$, and $D_{84}$ found from generating a new distribution is minimized compared against the actual values attributed to the segment in step 6.
9. Generate a distribution using the solved parameters, remove fractions finer than sand and coarser than the predicted $D_{max}$ 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.

Calculating $D_{16}$, $D_{50}$, and $D_{84}$ for measurements and quantifying uncertainty

For each measurement, the data is bootstrapped with replacement 100 times. For each bootstrap, values are selected 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_{16}$, $D_{50}$, and $D_{84}$ 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 contributing drainage area is recorded (from the drainage network attributes).

D max

for each measurement, the $D_{max}$ (i.e., the largest grain size measured) is also recorded in order to parameterize a function for estimating $D_{max}$ at each reach. If there is more than one measurement, a linear equation relating reach slope to $D_{max}$ is parameterized,

$D_{max} =\alpha S + \beta$

which can then be used at each segment. If there is a single measurement, The $D_{max}$ is divided by the $D_{84}$ value to record a factor by which $D_{84}$ can be multiplied to estimate $D_{max}$,

$f = \frac{D_{max}}{D_{84}}$

Critical Depth

For each measurement reach, after the $D_{16}$, $D_{50}$, and $D_{84}$ are calculated the depth required to mobilize each of these sizes (critical depth) is estimated using the Gilbert (2022) method for calculating critical Shields stress, along with the Shields function. Critical Shields stress for a given grain size is calculated as

$\tau_{c,i}^* = a\left ( \frac{D_{i}}{D_{50}} \right )^{-0.68}$

where the coefficient a varies as a function of relative bed roughness,

$a = 0.025$ if $\frac{D_{84}}{D_{50}} \leq 2$

$a = 0.087\ln \left ( \frac{D_{84}}{D_{50}} \right )^{-0.034}$ if $2 < \frac{D_{84}}{D_{50}} < 3.5$

$a = 0.073$ if $\frac{D_{84}}{D_{50}} \geq 3.5$

The critical Shields number is given by the function

$\tau^* = \frac{\tau}{\left ( \rho_{s} - \rho \right)gD}$

where $\tau$ is boundary shear stress, typically estimated using the depth-slope product (which assumes steady uniform flow). $g$ is acceleration due to gravity, $\rho_{s}$ is the density of sediment, and $\rho$ is the density of water.

$\tau = \rho g h S$

where $h$ is average flow depth and $S$ is reach slope. $\tau$ can be solved for by multiplying the denominator of the Shields function by $\tau^*$, and critical depth can then be estimated by rearranging the depth slope product,

$h_{c} = \frac{\tau_{c}}{\rho g S}$

Hydraulic Geometry Coefficients

In order to extrapolate critical depth values to each drainage network segment, coefficients for hydraulic geometry coefficients are parameterized. Hydraulic geometry equations scale channel geometry with contributing drainage area (Leopold & Maddock, 1953). Typically an exponent of 0.4 is used for depth:

$h = \alpha A^{0.4}$

Where $A$ is contributing drainage area. If there is a single grain size measurement, $\alpha$ is simply solved for that reach and used for every other drainage network segment. If there are multiple measurements, a function is parameterized to vary $\alpha$ with $S$, under the assumption that at steeper slopes, flow accelerates and depth is thus relatively more shallow for a given $A$. This function is a negative power function (i.e., increases in slope result in a non-linear decrease in depth)

$\alpha = a S^b$

Estimating Grain Size

After the function for calculating the coefficient in the depth hydraulic geometry formula is developed, critical depth for $D_{16}$, $D_{50}$, and $D_{84}$ is estimated for each reach, based on the reach's contributing drainage area. Shear stress at that depth is calculated using the depth slope product with the critical depth, and then the Shields formula is rearranged to solve for $D$.

$D_{i} = \frac{\tau_{c, i}}{\left ( \rho_{s} - \rho \right)g\tau_{c, i}^*}$

where the subscript, i represents the grain size percentile (i.e., 16, 50, 84).

The relationship developed for $D_{max}$ (i.e., predicting it as a function of slope or a factor of $D_{84}$ is then used to estimate $D_{max}$ at each reach.

Estimating Full Grain Size Distributions

For each grain size measurement, the data is converted by taking the square root of all measurements, then tested for normality. Typically, in coarse bedded, mountain streams grain size distributions are positively skewed, with a few high values from large boulders causing the skew (). The mean and standard deviation of the distributions is then computed. If there is more than one measurement, the mean of means, and standard deviations is taken.

For each network segment, a distribution is generated using the mean and standard deviation parameters from the measurement data. In an optimization process, those parameters are then varied until the error between $D_{16}$, $D_{50}$, and $D_{84}$ from a distribution generated from the new parameters and the values of $D_{16}$, $D_{50}$, and $D_{84}$ predicted for the segment previously is minimized. These optimized parameters are used to generate a final distribution, which is transformed back (by squaring it). Values finer than sand are removed, since these finer fractions readily move through mountain rivers and generally don't appear in the bed except in specific depositional zones. Additionally, any values coarser than the predicted $D_{max}$ for the reach are removed from the distribution. From the remaining values, the fraction within each half-phi size interval is then calculated and the values stored in a .json file where the format of the data is:

{

}

References

Leopold, L. B., & Maddock, T. (1953). The hydraulic geometry of stream channels and some physiographic implications (Vol. 252). US Government Printing Office.

Wolman, M. G. (1954). A method of sampling coarse river‐bed material. EOS, Transactions American Geophysical Union, 35(6), 951-956.