The Design of a Site-Calibrated Parker–Klingeman Gravel Transport Model

The use of site-calibrated models for predicting bedload transport in gravel-bed rivers remains relatively rare, despite advances in methodology and computing technology, and its notable advantages in terms of predictive accuracy. This article presents a new algorithm for site calibration of the Parker–Klingeman (1982) model, along with a detailed discussion of considerations that influence model selection and calibration methodology. New visualization techniques are explored to demystify the calibration process, using three examples with progressively more challenging calibration conditions. The new method is particularly well suited to streams with high sediment loads, or cases where extrapolation of transport function estimates is necessary.


Introduction
Accurately estimating bedload transport in gravel-bed rivers is difficult in practice. Transport equations are extremely sensitive, such that small errors in estimating hydraulic shear stress available to mobilize particles can lead to large errors in the predicted transport rate [1]. Accurately determining the characteristic grain size or grain size distribution of the active portion of the riverbed is also problematic [2]. Errors become amplified since grain size has a large effect on the predicted bedload transport volume, and since many models use riverbed grain size distribution to estimate the size distribution of the bedload. Finally, the thresholds for initiation of particle motion exhibit dependence on the way riverbed particles are either exposed to the flow or shielded by neighboring particles [3], which is, in turn, dependent on the distribution of grain sizes, and is thus also uncertain.
Early efforts to produce bedload models for practical applications relied on simple functions of stream power (e.g., [4]) or bed shear stress (e.g., [5]) in excess of a critical threshold for particle mobility. This critical shear stress, or critical Shields stress in non-dimensional form [6], was often assumed to be a constant, even though there were various interpretations of what that value should be, as summarized in the thorough review by Buffington and Montgomery [7]. Model progress advanced significantly when it was realized that the structure or arrangement of grains on the riverbed develops in response to establishment of a quasi-equilibrium between the bedload in motion and the riverbed material [3]. Thus, thresholds for grain mobility are not well represented by a single number, but depend on both the riverbed grain size distribution and the degree of surface coarsening [8,9].
In the decades since these fundamental developments, there have been efforts to improve on estimates of key model parameters. Andrews [8] predicted critical Shields stress for initiation of particle motion from the ratio of grain size to median grain size of the subsurface. Parker, Klingeman on a physically tenable bedload transport function, fewer samples are needed to achieve a given improvement in accuracy than would be the case if one were to fit the bedload sample data to a curve whose form is not known a priori, as is true with a purely statistical approach such as linear regression [19]. Furthermore, a site-calibrated model can be extrapolated more confidently than a statistical regression line, since it is based on physical principles that presumably extend beyond the range of available on-site data. Because the exact sources of error in estimating shear stress, initiation threshold of particle motion or grain-size distribution matter less when one uses a site-calibrated model, the technique can be used in streams that differ morphologically from the streams originally used to develop the transport function, and for a wider range of sediment sizes, including a sand component in mixed sand-gravel streambeds [20].
Our purpose here is to present a new site-calibrated model procedure, which will be shown to work well in situations that are hydraulically complex enough to render uncalibrated models impractical, but where only a small number of bedload samples exist. Advances also include the use of graphical techniques that facilitate an informed perspective of what the various calibration procedures are actually doing. Use of the model will be illustrated with three examples from study sites representing different degrees of modeling difficulty, as determined by bedload-data richness and hydraulic complexity: Oak Creek, Oregon, USA (a data-rich site, from which the Parker-Klingeman model [3,10] was originally developed), Paradise Creek, Oregon, USA (a site with a moderate amount of high-quality bedload data), and, finally, South Fork Thornton Creek in Seattle, Washington, USA (a site with a smaller amount of bedload data, which typifies many practical, but difficult, bedload modeling applications).
Finally, model calibration procedures greatly affect the applicability to particular modeling goals. These features include whether to optimize the fit of the model to bedload data using loads versus logarithms of loads, whether to include an adjustment to eliminate bias by matching total modeled load with total sample load, and the choice of representing the bedload by a single grain size, two grain size fractions or multiple grain sizes. A practical discussion of the implications of these decisions will be presented.

Materials and Methods
Our modeling procedure is based on the formulation of Parker and Klingeman [3], but site-calibrated differently than described in earlier literature [19][20][21].

Description of the Parker-Klingeman Model
The  model [3] is a multi-grain-size model, meaning that it predicts sediment transport for each of a series of grain size classes. It does this by virtue of a single bedload transport function, used in combination with a simple hiding factor equation. The hiding factor equation represents the effect of the coarsened surface layer (called pavement [3]) on the mobility thresholds for each grain size class.
More precisely, the hiding factor equation is where τ * ri = reference Shields stress (see below) for particles of size class i; D i = geometric mean particle diameter for size class i; D 50 = median particle diameter for the bed material layer used in the model (either surface or the underlying subsurface); τ * r50 = reference Shields stress associated with D 50 ; and γ = hiding factor exponent. The Shields stress for particle size class i is defined as Water 2017, 9, 441 4 of 16 in which τ = time-averaged bed shear stress, computed at a given point on the streambed as in Equation (3) and reflects local hydraulic conditions: Here, ρ = mass density of water; s g = specific gravity of the sediment particles; g = acceleration of gravity; d = local water depth; and S = energy slope.
The reference Shields stress is analogous to a critical Shields stress for initiation of particle motion, but is operationally defined as the value of Shields stress corresponding to a "small but measureable" sediment load [3], specifically, a dimensionless bedload rate of W * i = W * re f = 0.002, where In Equation (4), q Bi is the volumetric bedload transport rate per unit channel width, and f i is the fraction of the representative streambed sediment layer (either the pavement or subpavement) in grain size class i. As can be seen from the form of Equation (1), the hiding factor has the effect of reducing the reference Shields stress for particles larger than D 50 and increasing it for smaller particles. This tends to "equalize" the mobility [3,9] such that all grain sizes become mobile over a smaller range of discharges than would occur if there were no hiding factor effect.
The definition of the Shields stress ratio, Φ i , a ratio of Shields stress to reference Shields stress, allows the transport function for each individual grain size (see Figure 1) to fall onto the same curve: in which the constant β is initially set to 5600 (the value used by Parker and Klingeman [3]) and

Description of Model Calibration Procedure
Several coefficients in the above equations are potentially available for use in model calibration, including the reference Shields stress, τ * r50 , the hiding factor exponent, γ, and the transport function coefficient, β. In the calibration algorithm presented here, we have chosen to follow a procedure analogous to the original Parker-Klingeman (1982) model development [3] for determination of τ * r50 and γ, followed by an additional step for optimizing the value of β.
First, we compute the dimensionless bedload W * i and Shield stress τ * i for each grain size class of each bedload sample, using a cross section average shear stress corresponding to the sample conditions ( Figure 1). We then fit individual regression lines, by grain size class, through these bedload data. Since neither of these two variables can be assumed to be error-free, we used orthogonal regression [23] instead of ordinary linear regression for this step. The intersection of each of these lines with the reference value of W * i = 0.002 yields a reference Shields stress for the given size class, τ * ri . These τ * ri values are next plotted against a dimensionless grain size ratio, D i /D 50 , from which a corresponding regression line yields the hiding factor exponent, γ, and the D 50 reference Shields stress, τ * r50 (Figure 2). Since the bedload samples are integrated totals for the whole cross section, we used hydraulic radius (cross sectional area divided by wetted perimeter) instead of local depth to compute shear stress in Equation (3), from which Shields stress was calculated (Equation (2)) for each point in Figure 1. Next, the bedload is computed for each grain size at each point on the channel cross section for all samples from the foregoing model equations using the values of γ and τ * r50 determined previously, and using local depth to compute shear stress. These values are summed by sample, and a total for all samples together is computed. Finally, the constant β in the transport function (Equation (6)) is optimized such that the sum of total measured bedload, from all samples, equals the total computed bedload, using a simple ratio of total measured load to load computed initially with β set to 5600, as a multiplier.
Water 2017, 9, 441 5 of 16 equals the total computed bedload, using a simple ratio of total measured load to load computed initially with β set to 5600, as a multiplier.

Figure 1.
Size-specific transport rates. Log-log (power) regression lines are fit to the data for each grain size class to determine the size-specific reference Shields stress, * , the value of * where the fit line crosses * = 0.002. Regressions are computed to minimize sum of squared distances between the points and the regression line (orthogonal regression), except for the 32 mm line, which had too few points to accurately determine its slope. It was assigned a slope equal to the average of the 16 mm and 8 mm slopes, and was constrained to pass through the centroid of the data.  Size-specific transport rates. Log-log (power) regression lines are fit to the data for each grain size class to determine the size-specific reference Shields stress, τ * ri , the value of τ * i where the fit line crosses W * i = 0.002. Regressions are computed to minimize sum of squared distances between the points and the regression line (orthogonal regression), except for the 32 mm line, which had too few points to accurately determine its slope. It was assigned a slope equal to the average of the 16 mm and 8 mm slopes, and was constrained to pass through the centroid of the data.
Water 2017, 9, 441 5 of 16 equals the total computed bedload, using a simple ratio of total measured load to load computed initially with β set to 5600, as a multiplier.   Hiding factor equation. Reference Shields stress, τ * ri , for each grain size, determined from the intersection of the regression lines shown in Figure 1 with W * i = 0.002, are plotted against the grain size ratio D i /D 50 to obtain the hiding factor relationship. Here, D 50 refers to the subsurface layer. τ * r50 is found to be 0.155, and the exponent is 0.902. The first part of the procedure achieves minimum average sum of squared error (SSE) in dimensionless bedload by sample and grain size, while the second part, optimizing β, achieves zero bias, defined according to the calibration sample volume. The model was implemented using code written in the Python language, and was designed for easy input and output from Excel spreadsheets.

Application of Model Calibration to Three Study Sites
We tested the model calibration procedure using bedload and hydraulic data collected at three study sites representing a range of practical scenarios for bedload sediment transport modeling. Characteristics of the three sites are summarized in Table 1. Channel cross sections used in model computations are shown in Figure 3. Creek is a stream draining the east slope of the Coast Range near Corvallis, Oregon, USA. It has a drainage area of about 670 ha (2.6 mile 2 ), a slope in the study reach of about 0.014, a surface and subsurface D 50 of about 54 and 20 mm, respectively [24]. Peak flow hydrology is dominated by winter rain storms. Bedload samples were collected using a vortex sediment ejector extending the entire width of the channel [25]. Of 66 samples collected during the winter of 1971, 22 taken at discharges greater than 1 m 3 /s were used to develop the original Parker-Klingeman model [3,10,24]. These were re-analyzed for the current study. Based on the channel description in Milhous [24], the cross section shown in Figure 3c was used for modeling purposes.
Oak Creek, arguably, represents a "best case scenario" for sediment modeling. The dataset is large, and of high quality due to use of a fixed, full-channel-spanning sampler. The sediment load in Oak Creek was relatively stable, meaning that the relationship between bedload and discharge did not change appreciably over time. In fact, contemporary (2015-2016) measurements of bedload mass flux per unit width have been found to be consistent with those taken in 1971 [26], which suggests that bedload is strongly determined by transport conditions rather than exhibiting source dependence, and that reach hydraulic conditions and riverbed composition have been remarkably stable. The only minor complication was that the vortex sampler affected the local hydraulics to some degree, such that the average bed shear stress was not a uniformly increasing function of discharge. Moreover, use of a single water-surface-elevation to discharge relationship (rating curve) was not possible, and a separate equation was developed for the sequence of six samples beginning with No. 15 [24].  Creek is a stream in the headwaters of the Klamath Basin of South Central Oregon, USA, with drainage area of 6480 ha (25.0 mile 2 ). The watershed lies entirely within a National Forest, has a relatively stable sediment load over time, and snowmelt-driven peak flows that create relatively stable hydraulic conditions during bedload movement. Slope at the study reach is 0.0027, and surface and subsurface D50 are 39 and 18 mm, respectively [20]. Peak flows may occur during winter rain-onsnow events, but a consistent spring snowmelt pattern drives long-duration sediment transport. As part of a study to determine channel maintenance flows, 11 bedload measurements were taken during the snowmelt runoff season (April-May) of 1996, using a Helley-Smith sampler with an orifice of 76.2 mm. These same data were analyzed by Bakke et al. [20].
Paradise Creek represents a site with an intermediate amount of available bedload data. Ease of model calibration and certainty of bedload predictions are less than those encountered with the large Oak Creek dataset, but sufficient to serve as an accurate yet practical example of the calibrated approach. . The watershed lies entirely within a National Forest, has a relatively stable sediment load over time, and snowmelt-driven peak flows that create relatively stable hydraulic conditions during bedload movement. Slope at the study reach is 0.0027, and surface and subsurface D 50 are 39 and 18 mm, respectively [20]. Peak flows may occur during winter rain-on-snow events, but a consistent spring snowmelt pattern drives long-duration sediment transport. As part of a study to determine channel maintenance flows, 11 bedload measurements were taken during the snowmelt runoff season (April-May) of 1996, using a Helley-Smith sampler with an orifice of 76.2 mm. These same data were analyzed by Bakke et al. [20].
Paradise Creek represents a site with an intermediate amount of available bedload data. Ease of model calibration and certainty of bedload predictions are less than those encountered with the large Oak Creek dataset, but sufficient to serve as an accurate yet practical example of the calibrated approach. urbanized (with about 49% impervious surface), the creek flows through narrow, forested ravines cut into glacial drift of Pleistocene age, which serve as sources of gravel and sand bedload. The system responds rapidly to winter rainstorms, which produce the peak flows of record. The South Fork has a watershed area of 910 ha (3.5 mile 2 ), a surface and subsurface D 50 of 39 and 16 mm, respectively, and a slope at the study site of 0.012.
As part of design planning for a habitat restoration project, a calibrated bedload model was developed for the study reach near the confluence with the North Fork of Thornton Creek, utilizing seven bedload samples that were collected between November 2009 and March 2012 [27]. Some of the samples were collected with a 76.2 mm Helley-Smith bedload sampler and some with the larger Elwha sampler (101.6 by 203.2 mm opening). A sediment detention pond on the mainstem of Thornton Creek just downstream from the confluence provided an empirical check on long-term sediment load estimates.
South Fork Thornton Creek represents a difficult modeling exercise, but one typical of many sites where bedload estimates are needed. The mainstem of Thornton Creek, based on detention pond dredging records, has an annual sediment load of 27 to 31 m 3 /km 2 ; most of this sediment originates in the South Fork. For comparison, Thornton Creek has more than double the pre-development annual sediment load of 12 m 3 /km 2 for the nearby Issaquah Creek watershed [28]. This relatively high sediment load, combined with rapidly changing water levels during peak flows, amplifies the potential for error in estimation of cross section, slope, bed material and hydraulic characteristics. For these reasons, application of a theoretical (uncalibrated) model at this site would be problematic.

Results
Results of model calibration for the three study sites are given in Table 1, which also summarizes site characteristics and modeling challenges. Total bedload is shown in Figure 4, using cross-section average shear stress as the independent variable. Since the use of bedload sediment rating curves is common and intuitive, we also plot the model calibrations results with water discharge as the independent variable in Figure 4.
Comparison of the bedload plotted against average shear stress versus water discharge for each site reveals less scatter when shear stress is used as a predictor of bedload. This is particularly evident at Oak Creek, where six of the samples were taken during hydraulic conditions that necessitated a different water surface elevation-to discharge rating curve. Inflexions in the lines representing model-predicted bedload versus shear stress for Paradise Creek and South Fork Thornton Creek are due the fact that these lines were determined using equally-spaced intervals of water discharge, back-calculating the water surface elevation from the discharge rating curve, and then computing an average hydraulic radius and corresponding average bed shear stress. Bedload, however, was computed at each of 25 or more locations on the channel cross section using local depth to compute shear stress, and then summed to give total bedload. Since bedload increases exponentially with depth, inflexions in the cross section shape, and inflexions due to the use of a compound rating curve (necessary with overbank flow), become amplified, and total computed bedload is influenced more greatly by zones of maximum depth than by average shear stress. Figure 5 shows how the sum of squared error (SSE), determined from the calibrated model and calculated as the square of [log(measured load) − log(computed load)], varies with the choice of reference Shields stress (τ * r50 ) and hiding factor exponent (γ). Each grain size for each sample is included in this sum. Two observations require further explanation.
First, the use of log units equalizes the influence of large and small samples, but results in the "spiky" appearance of these graphs. The spikes in the figures represent locations in the field of τ * r50 versus γ where a grain size is at the threshold of incipient motion, according to the model transport function. Since the logarithm of zero is undefined, this grain size had to be excluded from the computation of squared error. Moving slightly in either the τ * r50 or γ axis direction, one reaches a location where a miniscule amount of bedload of that grain size is predicted. In logarithms, that bedload corresponds to a large negative number, resulting in the abrupt appearance of a "large" spike in SSE.
Secondly, the optimum combination of τ * r50 and γ, as determined from the method described earlier using reach-averaged shear stress and total bedload for the whole cross section, does not necessarily fall on the apparent lowest SSE point in Figure 5. This is due to determination of optimum τ * r50 and γ from average shear stress, while the model is computing bedload (and thus SSE) from locally-determined shear stress along the cross section.   First, the use of log units equalizes the influence of large and small samples, but results in the "spiky" appearance of these graphs. The spikes in the figures represent locations in the field of * versus γ where a grain size is at the threshold of incipient motion, according to the model transport function. Since the logarithm of zero is undefined, this grain size had to be excluded from the computation of squared error. Moving slightly in either the * or γ axis direction, one reaches a location where a miniscule amount of bedload of that grain size is predicted. In logarithms, that bedload corresponds to a large negative number, resulting in the abrupt appearance of a "large" spike in SSE. In the case of the most difficult example presented herein, South Fork Thornton Creek, the new Parker-Klingeman model calibration algorithm outperformed several similar approaches attempted ( Figure 6). For example, the algorithm of Bakke, et al. [20] is calibrated by adjusting the hiding factor exponent to obtain minimum logSSE while maximizing grain movement prediction accuracy, and then adjusting the reference Shields stress (τ * r50 ) to obtain zero bias. This algorithm could be calibrated easily, but the resulting transport curve has a very large slope, and predicts unreasonably large transport rates when extrapolated to larger peaks flows, unlike the new algorithm described above. Bedload Assessment for Gravel-bed Streams (BAGS), a publicly available spreadsheet-based modeling package [21], which calibrates the Parker-Klingeman model by minimizing arithmetic SSE, would not converge to a calibrated solution for this site, probably due to the high sediment loads measured. Theoretical models, including the (uncalibrated) version of Parker-Klingeman, all yielded unreasonably high bedload transport predictions, when juxtaposed against the total annual load estimates from dredging records [27].  Table 1.
Secondly, the optimum combination of * and γ, as determined from the method described earlier using reach-averaged shear stress and total bedload for the whole cross section, does not necessarily fall on the apparent lowest SSE point in Figure 5. This is due to determination of optimum * and γ from average shear stress, while the model is computing bedload (and thus SSE) from locally-determined shear stress along the cross section.
In the case of the most difficult example presented herein, South Fork Thornton Creek, the new Parker-Klingeman model calibration algorithm outperformed several similar approaches attempted ( Figure 6). For example, the algorithm of Bakke, et al. [20] is calibrated by adjusting the hiding factor exponent to obtain minimum logSSE while maximizing grain movement prediction accuracy, and then adjusting the reference Shields stress ( * ) to obtain zero bias. This algorithm could be calibrated easily, but the resulting transport curve has a very large slope, and predicts unreasonably large transport rates when extrapolated to larger peaks flows, unlike the new algorithm described above. Bedload Assessment for Gravel-bed Streams (BAGS), a publicly available spreadsheet-based modeling package [21], which calibrates the Parker-Klingeman model by minimizing arithmetic SSE, would not converge to a calibrated solution for this site, probably due to the high sediment loads measured. Theoretical models, including the (uncalibrated) version of Parker-Klingeman, all yielded unreasonably high bedload transport predictions, when juxtaposed against the total annual load estimates from dredging records [27].  Table 1.

Discussion
Several modeling considerations are illustrated by these examples. Although these issues are particular to the site-calibrated modeling approach, they are inherent in the development of published theoretical models as well. In the site-calibrated approach, these considerations are explicitly addressed, which ultimately improves confidence in the results.

Optimization Using Arithmetic Loads versus Logarithms of Loads
First, development of a site-calibrated model requires choice of criteria for fitting the model to  [3] without calibration (red line). Paradise Creek, shown on the left, has relatively low shear stress during bedload transporting events compared to the South Fork Thornton Creek, on the right, resulting in a striking difference between uncalibrated versus calibrated model performance.

Discussion
Several modeling considerations are illustrated by these examples. Although these issues are particular to the site-calibrated modeling approach, they are inherent in the development of published theoretical models as well. In the site-calibrated approach, these considerations are explicitly addressed, which ultimately improves confidence in the results.

Optimization Using Arithmetic Loads versus Logarithms of Loads
First, development of a site-calibrated model requires choice of criteria for fitting the model to the data. One way to do this is to use a minimum average of squared errors (SSE) approach, where errors are defined as the difference between the computed and measured loads. This is the approach used in BAGS [21]. A disadvantage of using arithmetic loads to compute SSE is that the model calibration will be dominated by the larger samples. This can result in unrealistically high slopes for the load versus shear stress relationship. This can be a problem if the model is intended to be extrapolated for computation of load at large discharges, as would be the case in use for an effective discharge computation, or for computing average annual load [29]. Furthermore, for some modeling objectives, such as determining threshold discharges for grain mobilization, or predicting sediment loads or size composition under conditions of marginal transport, the smaller samples should be weighted equally to the larger samples in the modeling optimization process.
One simple way to do this, which was the procedure followed here (as well as by Bakke et al., [20]), is to compute SSE from the logarithms of loads, rather than the actual loads. However, since the logarithm of zero does not exist, combinations of model parameters that yield zero load for some of the grain samples will need to be excluded from the computation. Moreover, combinations of parameters that produce tiny amounts of predicted sediment transport will result in "spikes" in the two-dimensional surface representing SSE, since the logarithm of this tiny amount will be large and negative.
The exclusion of zero-load grain samples creates a dilemma, in that a "low" SSE region in the Figure 5 can develop due to exclusion of grain samples, rather than best model fit. Thus, another criterion, such as optimum proportion of correctly predicted movement or non-movement of sample grain sizes, or zero overall bias, needs to be added to in order to interpret the log-load SSE diagram for model best fit. As an example, Figure 7 shows grain movement prediction accuracy for South Fork Thornton Creek.

Rationale for a Two-Stage Optimization
In our calibration procedure, optimization is accomplished in two steps. First, in development of the hiding factor constants τ * r50 and γ, SSE in dimensionless bedload, W * i , and Shields stress, τ * ri , are minimized to develop the regression lines used to predict τ * ri . This is done in log units, but since only the (finite) bedload samples are used, there is no zero exclusion or spiking issue. Second, the transport function constant is adjusted to achieve equality between total sampled and computed bedload (zero bias). Alternatively, this could have been done with a second stage of SSE minimization, using W * i versus Φ i , as Parker and Klingeman [3] originally did to fit the transport function curve to their data. Use of the zero-bias approach, however, is simpler. It also effectively adapts the original transport function, which was based on reach-averaged shear stress, to the quasi-two-dimensional computation being used here, which computes bedload at each point on the cross section.
for model best fit. As an example, Figure 7 shows grain movement prediction accuracy for South Fork Thornton Creek. Creek. Prediction error is the count of all grain sizes, in all samples, where the model predicted sediment movement when none was measured, or where sediment was measured but none was predicted. Smaller numbers signify greater accuracy in predicting grain movement. There were a total of 52 measurements by grain size in the seven available bedload samples.

Rationale for a Two-Stage Optimization
In our calibration procedure, optimization is accomplished in two steps. First, in development of the hiding factor constants * and γ, SSE in dimensionless bedload, * , and Shields stress, * , are minimized to develop the regression lines used to predict * . This is done in log units, but since only the (finite) bedload samples are used, there is no zero exclusion or spiking issue. Second, the transport function constant is adjusted to achieve equality between total sampled and computed bedload (zero bias). Alternatively, this could have been done with a second stage of SSE minimization, using * versus , as Parker and Klingeman [3] originally did to fit the transport function curve to their data. Use of the zero-bias approach, however, is simpler. It also effectively adapts the original transport function, which was based on reach-averaged shear stress, to the quasitwo-dimensional computation being used here, which computes bedload at each point on the cross section.

One-Dimensional versus Quasi-Two-Dimensional Model Approach
Typically, sediment transport models are made one-dimensional, meaning that the sample site is represented by a single average cross section and average shear stress. This approach works well for streams like Oak Creek, which has a flume-like, rectangular cross section, but not as well for Paradise Creek, South Fork Thornton Creek or other typical alluvial streams that have a definite Thalweg. Local shear stress can be quite different in the deeper Thalweg than the cross-section average, which has implications for model prediction accuracy. A simple way to account for this is to compute shear stress and bedload using local depth rather than average hydraulic radius, and to sum the incremental bedload values over the cross section. Since the bedload transport function is highly non-linear, the difference in computed bedload using this quasi-two-dimensional approach versus a single average shear stress can be striking, and leads to optimum values for the transport function coefficient, β, which are quite different from the literature. It also results in more accurate prediction of the largest grain sizes moved for a given total bedload. Creek. Prediction error is the count of all grain sizes, in all samples, where the model predicted sediment movement when none was measured, or where sediment was measured but none was predicted. Smaller numbers signify greater accuracy in predicting grain movement. There were a total of 52 measurements by grain size in the seven available bedload samples.

One-Dimensional versus Quasi-Two-Dimensional Model Approach
Typically, sediment transport models are made one-dimensional, meaning that the sample site is represented by a single average cross section and average shear stress. This approach works well for streams like Oak Creek, which has a flume-like, rectangular cross section, but not as well for Paradise Creek, South Fork Thornton Creek or other typical alluvial streams that have a definite Thalweg. Local shear stress can be quite different in the deeper Thalweg than the cross-section average, which has implications for model prediction accuracy. A simple way to account for this is to compute shear stress and bedload using local depth rather than average hydraulic radius, and to sum the incremental bedload values over the cross section. Since the bedload transport function is highly non-linear, the difference in computed bedload using this quasi-two-dimensional approach versus a single average shear stress can be striking, and leads to optimum values for the transport function coefficient, β, which are quite different from the literature. It also results in more accurate prediction of the largest grain sizes moved for a given total bedload.
In electing to use this approach, another issue arises, however. Since the bedload samples are composites for the whole cross section, there is insufficient information to compute W * i values for individual points on the cross section. Only an average W * i can be computed from the data, and this requires a corresponding single (e.g., average) value of Shields stress τ * i for determination of the hiding factor coefficients. However, when the model computes bedload according to the transport function (Equation (6)), a local value of Shields stress is used. The result of this difference is that the bedload transport rate, q Bi , corresponding to W * i = W * re f = 0.002, is different when computed from local depth than what its corresponding value would be when based on cross-section average shear stress in the hiding factor coefficient analysis. The second stage of the calibration process, adjusting the transport function coefficient β, eliminates this difference. Although adjustment of β could be done by another round of minimization of SSE, shifting the transport function vertically downward to pass through either the centroid of the data or to equalize the sum of total computed and measured loads is far simpler as it eliminates the complexity associated with the zero exclusion issue without affecting the fit of the curve to the data in a substantial way. Moreover, we opted to use the zero bias criterion as a way of insuring that cross sections with deep Thalwegs would not yield models that overpredict bedload at larger discharges. If predictions at marginal transport were the main objective, then an approach of adjusting β such that the transport function intersects the centroid of the data would be an appropriate choice.

Advantages of a Multi-Grain-Size Model
More generally, another issue that arises in transport modeling is which model to use, and whether to use a model that predicts total sediment volumes only, or the grain-size distribution of the bedload sediment. In regard to the data-calibrated approach described herein, two considerations are paramount.
First, when only a few bedload samples are available, as was the case here with South Fork Thornton Creek, practitioners need tools for assessing data quality. One of the advantages of a multi-grain-size model is that the grain size distribution of the data contains useful information about its appropriateness for model calibration. If the site fits the assumption that the bedload approaches an equilibrium with the streambed material, which is implicit in all of the physically based models, then the data should produce a family of nearly-parallel curves such as displayed in Figure 1, and the slopes of these curves should be close to that expected for a bedload transport function, which is 4.5 in the case of the Parker-Klingeman model. The information contained in the relationship between grain sizes thus becomes part of the calibration process, which effectively expands a single bedload measurement into a sub-set of measurements, one for each grain size. This allows the practitioner to spot departures from equilibrium and, conversely, to spot irregularities in data that might suggest poor sample quality if equilibrium is expected due to other lines of evidence. The patterns visible in bedload data stratified by grain size help the modeler to diagnose whether a sample "outlier" represents a measure of the variability found in nature or, conversely, represents a sample that is deficient in some regard, and should justifiably be excluded from the calibration. An example of this departure is shown in Figure 8, for the North Fork of Thornton Creek, which was a sampling site in close proximity to the confluence with the South Fork presented above. This diagnostic power is not available when using a model that predicts only a single, average sediment volume or two components (gravel, sand), as opposed to a series of multiple grain sizes. Samples marked with open circles are not consistent with the rest of the data due to differing site hydraulics, and were eliminated from the calibration process. This example is provided to show how the multigrain-size modeling approach can be used to check data quality and sediment equilibrium assumptions.
Finally, although most any model can, in principle, be calibrated with bedload data, the most state-of-the-art gravel transport models incorporate some form of a "hiding factor" to account for the way that the structure of the streambed causes particle mobility by grain size to adjust over what it would be in a streambed of uniform-sized particles, and thereby to achieve equilibrium between the streambed material and the bedload in transport [3]. Without this hiding factor, which effectively increases the apparent mobility of larger particles and reduces that of smaller particles, a model will invariably over-predict transport of small grain sizes and under-predict the large sizes, and typically will over-predict the total sediment load [22]. Moreover, a site-calibrated model should incorporate a transport function whose form derives from basic physical principles, ensuring consistency with Samples marked with open circles are not consistent with the rest of the data due to differing site hydraulics, and were eliminated from the calibration process. This example is provided to show how the multi-grain-size modeling approach can be used to check data quality and sediment equilibrium assumptions.
Finally, although most any model can, in principle, be calibrated with bedload data, the most state-of-the-art gravel transport models incorporate some form of a "hiding factor" to account for the way that the structure of the streambed causes particle mobility by grain size to adjust over what it would be in a streambed of uniform-sized particles, and thereby to achieve equilibrium between the streambed material and the bedload in transport [3]. Without this hiding factor, which effectively increases the apparent mobility of larger particles and reduces that of smaller particles, a model will invariably over-predict transport of small grain sizes and under-predict the large sizes, and typically will over-predict the total sediment load [22]. Moreover, a site-calibrated model should incorporate a transport function whose form derives from basic physical principles, ensuring consistency with the body of work under which the original model was derived.

Conclusions
We have shown that predicting bedload sediment flux using site-calibration of a bedload transport model can achieve high accuracy using only a small number of flux measurements. In particular, we developed a new algorithm for site calibration of the Parker-Klingeman (1982) model, which predicts bedload flux by particle size class. The approach is applied in three study sites representing a spectrum of hydraulic complexity and data availability, and is shown to be especially well suited to streams with high sediment loads, and to applications where extrapolation of the predicted sediment rating curve is required. Key elements underlying the success of this approach include: log-transformation of sediment loads to optimize the model for minimum error; accounting for lateral variation in shear stress; and the use of a multi-grain-size model, which maximizes the information value of each flux measurement, and by including a 'hiding function', accounts for the dynamics of particle interactions that govern the stress-flux relationship. Finally, we provided a novel graphical depiction of the calibration process, to assist practitioners in understanding and applying the site calibration technique.