DEVELOPMENT OF ONE-DIMENSIONAL BED VARIATIONANALYSIS METHODS FOCUSED ON THE MECHANISM OF
SEDIMENT TRANSPORT IN STONY-BED RIVERS
Kengo Osada1 and Shoji Fukuoka2
1

Research Associate, Research and Development Initiative, Chuo University
1-13-27 Kasuga, Bunkyo-ku, Tokyo, 112-8551, Japan, e-mail: osada@tamacc.chuo-u.ac.jp
2
Professor, Research and Development Initiative, Chuo University
1-13-27 Kasuga, Bunkyo-ku, Tokyo, 112-8551, Japan, e-mail: sfuku@tamacc.chuo-u.ac.jp

ABSTRACTIn a stony-bed river with wide grain-size distribution, boulders that resist movement
during a flood act as strong fluid-resistance element. In addition, sand and gravel often remainin the bed irregularities formed by large stones. Therefore, a stony-bed river will remain in astatic equilibrium condition, with no sediment discharge, when the river bed surface grain size
distribution corresponding to tractive force of the flow are reached. However, conventionalriverbed variation analysis does not account for this essential mechanism in stony-bed rivers.
In this study, the authors developed a new method of one-dimensional (1-D) bed variation
analysis for stony-bed rivers with well-graded sediment by considering the riverbed’s stability
mechanism created by large stones. We verified the applicability of the model with data fromfield experimentation carried out in the Joganji River.
Keywords: stony-bed river, well-graded sediment distribution, field experiment,
bed variation analysis, boulder, static equilibrium

1.

INTRODUCTION

Much research has been conducted to elucidate the mechanism of bed variation instony-bed rivers—rivers in which the bed material is comprised of stone and gravel having alarge grain-size distribution (e.g., Thorne,C.R. et al. 1987). However, most results obtained
pertain to riverbeds comprised of gravel having particle sizes ranging from 0.2 to 7.5 cm;
factors such as the difficulties of measurement in swiftly moving flows have impeded
adequate research on the mechanism of bed variation in stony-bed rivers containing largestones (i.e., >30 cm).
The authors (Kuroda and Fukuoka et al. 2005; Fukuoka et al. 2006) have previously
conducted large-scale field experimentation in the Joganji River to elucidate the mechanismof bed variation in stony-bed rivers. In a stony-bed river, bed scouring caused by flood flows
exposes boulders and other large-diameter bed material, the immovability of which acts as alarge resistance element against flood flows. Furthermore, the shielding effect of these large,
immovable stones causes small- and medium-diameter gravel to collect in the wake of theseobstacles. Consequently, the formation in a stony-bed river of a bed shape and bed material
size distribution suited to the flow will achieve static equilibrium in which neither stone nor
gravel are transported. This is an important mechanism of bed variation in stony-bed rivers.
Conventional techniques for analyzing bed variation in stony-bed rivers in generalconsist of bed load equation with mixtures and continuity equations for sediment and grain

sizes (Hirano, 1971). The equations of Ashida and Michiue (1972) are widely used to
calculate bed load equation for sediment mixtures. To calculate critical tractive force forsediment mixtures in those equations, the equation of Egiazaroff (1965) as modified by
Ashida and Michiue (1972) are often used. However, Fukuoka et al. (2006) have compared
their results with critical tractive force for sediment mixtures calculated with the modified
Egiazaroff equation and found that in an actual bed that has achieved static equilibrium, themodified Egiazaroff equation indicates that gravel would be transported. This showed that themodified Egiazaroff equation does not adequately account for the aforementioned shielding
effect of large stones and so is difficult to apply to stony-bed rivers.
In a stony-bed river, stones and gravels are moved violently during flooding, readilyleading to bank erosion and bed scouring near revetments. Countermeasures to these
problems include underpinning for existing revetments and the construction of new
revetments at locations of bank erosion, but such measures have been ineffective because of
the inability to predict the location and scale of scouring and erosion. In addition, dams andother lateral structure have blocked the downstream transport of large-diameter bed material,
while materials of specific sizes have been removed from riverbeds for use in construction,etc., thereby altering grain-size distributions. The loss of large-diameter materials has reduced
the ability of such riverbeds to resist flood flows, which could lead to even greater bed
scouring. To determine which bed material populations contribute to the formation of stablebeds, it is important to establish a highly accurate technique for analyzing bed variation instony-bed rivers.
The authors have established a one-dimensional bed variation analysis method that
focuses on the mechanism of stone and gravel transport. This technique takes an integrated
view of the temporal changes in water surface profile and bed variation and how thesequantities vary according to discharge. In particular, the determination of bed variationaccounts for the static equilibrium of bed materials, which is an important mechanism of bed
variation in stony-bed rivers. This technique was validated by applying it to large-scale field
experimentation in the Joganji River and comparing the resulting calculations to the observed
data. This paper also discusses the differences in approach between the authors’ technique andexisting techniques and enumerates the problems with the latter.

2.

A NEW METHOD FOR THE ONE-DIMENSIONAL ANALYSIS OF BED
VARIATION

This section describes the new technique that the authors have formulated for the
analysis of bed variation in stony-bed rivers—a method that takes into account the mechanismof stone and gravel transport. The underlying philosophy of this method is also discussed.
In a stony-bed river containing large stones, these stones are exposed in the course oferosion and thereby create resistance to the fluid field. In addition, the stones’ shielding effectcauses gravel to collect around the stones. Because large stones thus contribute greatly to bed
stability, methods for analysis must adequately account for this effect. Believing that theeffect of large stones should be incorporated as a form of bed irregularities, and rather thanassuming constant bed height in the computational grid as conventional bed variation analysis
techniques do, the authors devised a method of calculating the height of each size group ofmaterials on the riverbed surface. To incorporate this method, the authors created equations inwhich a transported-gravel layer (i.e., movable particle layer) is completely separated fromthe layer of nonmovable stones and gravel (i.e., the bed) (see Figure 1) and sediment
discharge and bed variation are calculated using quantities picked up from the bed anddeposited to the bed. The flow of calculations is shown below.

ε d 90

Figure 1 Concept of the new analysis method

Figure 2d 90 (m )Figure 3 Relationships between ε d 90 and d 90

Process of the new bed variation analysis

Figure 2 shows the flow of calculations in the authors’ bed variation analysis
technique. The flow field is calculated with equations (1) and (2). Although the bed resistanceterm generally used in the equation of motion is Manning’s roughness coefficient or someother representation of frictional resistance, in a river with stones and other large-diameterbed materials, the form resistance of these materials is predominant. Therefore the authors,
viewing riverbed materials of size d80 or greater as contributing to resistance (Fukuoka et al.,2006), determined resistance with equation (3), which uses d90, the intermediary between d80
and the maximum diameter. Nd90 is the number of particles contributing to bed resistance perunit area. As particles of diameter d80 or greater contribute to resistance, this number is 20%

of the unit area. The number of d90 in 20% of the unit area is calculated with the following
equation.

N d 90 =0 .2
α 2 d 902

(14)

where α2,α3 = 2-D and 3-D shape factors of the particles, CD = drag coefficient. Since actual
gravel particles have an irregular shape, not a spherical one, and since the particles are
assumed to be half-buried, a value of 1.0 was used for CD in equation (3), which is also basedon the values for submerged roughness drag coefficient observed by Uchida and Fukuoka etal. (2001). In equation (3), εd90 is the shielding coefficient of d90 and the product of εd90α2d902 is
the projectional area of d90 in the flow direction. Figure 3 shows the correlation between theexperimentally determined (Kuroda and Fukuoka et al., 2005) d90 and εd90. Based on this, theequation for calculating εd90 is as follows.
1.5

d 
ε d 90 = 0.35 90 (15)
 0.3 Also, uf90 is the velocity acting on d90 and is determined with a logarithmic velocity
distribution equation.

(a) Concept of the pick up particle(b) Calculation method of γ k ,k ′Figure 4 Concept of the pick up particle and calculation method of γ k ,k′The first step in bed variation analysis is to determine whether or not each particle sizegroup is picked up according to balance-of-moment equation (equation (4)). An overview ofthe process of pick up judgment is shown in Figure 4(a). For each particle whose pick up is tobe determined, a stationary particle (of diameter d80) is established below angle θk, andequation (4) is calculated with its point of contact as the fulcrum. Here, CL is the lift
coefficient, which was assigned a value of 0.2 based on the work of Fukuoka and Watanabe etal. (2005). In addition, αx is the distance in direction x from the contact point to the center ofthe separated particle, αz the distance in direction z from the contact point to the point where
fluid force acts, and uf the velocity acting on the particle, which are calculated according tologarithmic law using the average bed (discussed below) as the standard surface. The
subscript k indicates the particle diameter group. The point where fluid force acts was definedas the 0.8dk portion of a separating particle, as shown in Figure 4(a). The angle θk was
calculated as shown below using diameter d80.

θ k = 90.0 − 45.0 exp (Z B i , k − Z d 80 i ) 80 i 
d

(16)
dk 
where ZBi,k = height of diameter dk, Zd80i = height of d80, and subscript i indicates thecomputational grid number. This equation represents the relationship between d80 and theheight–diameter ratio of the picked up particles. For instance, when the diameter is less thand80 and height is low, then angle θk is greater than 45°, indicating pick up is less likely.

Equation (16) shows that the shielding effect of large-diameter particles reduces the pick up
rate of smaller-diameter particles and is also relevant to pick up time, which is discussed ingreater detail below. This is an important indicator of pick up, although because sufficientverification has not been conducted, additional study is necessary.
Next we explain the method of determining bed stability. Equation (5) represents theconditions for bed stability. The left side is the downstream-direction component of fluidmass, while Fx' is the resistance from those bed particles of size d80 or greater determined byequation (4) to resist transport. We assume that when all particles of size d80 or greater resist
transport, it means that all drag from the fluid field is absorbed by those particles, and fluid
force does not act on the surrounding particles. In the calculations, pick up analysis isconducted on particles of size d80 and d90; when both are found to resist pick up (i.e., whensizes d80 and greater are stable), it is assumed that all fluid drag is absorbed by particles ofsize d80 and greater, with no drag acting on the other particles, and so the quantity of pick up
from the bed is deemed to be zero.
After the process outlined up to this point, equation (6) is used to determine thequantity of pick up , i.e., the quantity of particles found to pick up from the bed. Here, VPi,k isthe pick up quantity per unit time and unit area, NPi,k is the number of particles of eachdiameter group in the bed, and TPi,k is the time required for pick up from the bed—using theequation of motion for particles, this is defined as the time required for a particle to betransported to the top of a stationary particle, as shown in Figure 4(a). However, even if pickup judgment shows particles of diameter k to be picked up, it is possible that not all particles
of diameter k on the bed surface will pick up. For instance, some such particles may beshielded by larger particles, or irregularities in the bed surface may prevent forces from acting
uniformly on all particles of that diameter. Consequently, rather than picking up all particles
of diameter k on the surface, the number of picking up particles must be controlled. In
equation (6), εPi,k is a coefficient that accounts for the shielding effect of particles larger thandk, with no shielding effect on particles larger than d80 (εPi,k = 1.0), while the shielding effecton particles smaller than d80 is calculated with the following equation.
k −1

ε P i , k = 1.0 − ∑ (γ k ,k ′ ⋅ Pi ,k ′ )(17)

k ′ =1

Pi,k is the bed surface ratio of particles of diameter dk, and γk,k′ is defined as follows in relation
to the relative heights of particles k and k′. As shown in (1) in Figure 4(b), the height of thebottom of particle k is compared to the height of the top of particle k′, and if the bottom of k ishigher, then the value of γk,k′ is zero, i.e., particle k′ has no shielding effect on particle k. When,as in (2) in Figure 4(b), the height of the bottom of particle k is compared to the center ofparticle k′, if the bottom of particle k is lower then the value of γk,k′ is 1. For particles between(1) and (2) the value of γk,k′ is determined linearly as in Figure 4(b). PPi is the ratio between
the fluid field force acting on the bed and the sum total of drag acting on the particle, and is
calculated as follows:

PP i = nρghi I e(18)
 N P i ,k2 2 ρCDα 2 d k u f ∑2
k =1 

Next, we explain how deposition rate on the bed is calculated. The amount ofdeposition on the bed is calculated with equation(3), in which VDi,k is bed deposition quantity
per unit time and unit area, PCi,k is the collision rate (i.e., rest rate), and VSi,k is the sedimentvolume within the computational grid, which is calculated with equation (8). In addition, Ai is
the computational grid area ( Ai = ∆x i × Bi ) and upi,k in equation (8) is particle velocity in thedownstream direction.

The rest ratio is an important variable in the equation for bed deposition. The author’s
method includes calculating 5-second saltation using the particle equation of motion, thendetermining the bed particle collision ratio, which is assumed to be the rest ratio (Figure 5).
An actual riverbed is comprised of a broad range of particle size populations and has
considerable irregularities, with particles conceivably being stopped by depressions andcollisions with large particles. However, the authors’ method does not yet account for this.
The method for determining the rest ratio must be improved through the detailed examinationof field data. Particle velocity (upi,k) is determined by dividing the calculated time (5 seconds)
into the particle’s downstream movement distance as determined by the saltation time
calculations. In addition, sediment discharge is calculated not with the sediment-transport
equation but in a way that accounts for the balance between pick up rate and deposition rate.

Figure 5

Concept of the saltation analysis

Next, let us see how the authors calculate the height of each particle size group(defined as the height of the particle tip), bed surface ratio of dk particles and average bedheight. The height of each particle size group moment by moment is calculated using eachsize group’s pick up quantity from the bed, deposition quantity to the bed, and the bed surface
ratio of dk particles with equation (10). As seen in equations (11) and (12), the bed surface
ratio of each particle size group is calculated using each size group’s pick up and deposition
quantity. Here, P0i,k is the subsurface particle diameter ratio. Average bed height is needed to
calculate the flow field and is calculated as shown in equation (13). The first term on the rightside of equation (13) is the average particle height in the computational grid, and the second
term on the right is radius of the average-diameter particle. Thus, average bed height isdefined as average particle height less the radius of the average-diameter particle.

3.

VALIDATION THROUGH APPLICATION TO 2004 FIELD EXPERIMENTS

3.1 Conditions of Calculations
The authors tested their bed variation analysis method by applying it to the results offield experiments conducted in the Joganji River in 2004. In the experiment, a steady flow
was released through a 170-meter-long linear waterway. Water level, discharge, bed height,
and bed surface particle size distribution were measured. Figure 6 shows the profile of thewaterway’s bed height, water level, and waterway width. A rapidly-flowing segment formedin the waterway’s middle, where width also narrowed. This narrowing occurred because prior
to the start of the experiment, bed materials that were larger than materials elsewhere in thewaterway were in existence there. Consequently, the initial surface and subsurface particle
size distribution used for this swiftly-flowing section in the calculations was larger than theslowly-flowing section. The grain size distribution was based on the dimensionless particle

Figure 6 Longitudinal distribution of water
elevation, bed elevation and waterway width

Figure 7 Initial condition of grain
size distributions

size distribution determined by Kuroda and Fukuoka et al. (2005). Seven particle sizes wereused in the calculations—35 cm, 27 cm, 20 cm, 12 cm, 8 cm, 5 cm, and 2.5 cm. The initial
bed height used in the calculations was that from the experiment (Figure 6), and the bed
variation calculations were run until bed height stabilized throughout the entire waterway.
3.2 Results of Calculation
Figure 8 compares the calculated and experimentally obtained water level and bed
height profiles over time; Figure 9, bed transport rate over time of each bed material size. In
the experiment, static equilibrium was achieved approximately 24 minutes after the flow ofwater began. In the calculations, bed material transport had nearly stopped—excepting in thedownstream section—after 30 minutes and within 40 minutes had completely stopped
throughout the waterway, achieving static equilibrium. Looking at the results for water leveland bed height after 40 minutes (Figure 8(c)), we see that the calculations accurately
reproduce the experimentally obtained water surface profile and bed height. Figure 8(c) also
shows the results of calculations performed with a conventional technique (described below)
simulating conditions after 5 hours of flow. In these conventional simulations, large-diameterbed particles were unmoved while small-diameter particles moved continuously, resulting in
extremely slow bed variation and a bed height that continued to drop even after 5 hours. In
short, results differed greatly from the experimentally observed values for both water leveland bed height. The data in Figure 9 for the transport of each bed material size indicate
extensive transport of 8-cm, 12-cm, and other cobble-class material, whereas transport of 35cm, 27-cm, and other large-sized rocks was limited. The low amount of transport of large
rocks is probably due to the low pick up amount and the slow transport velocity, while theshielding effect of the large-sized bed material—which was accounted for in thecalculations—reduced the pick up rate of 2.5-cm and other small-sized bed material. Figure
10, a comparison of experimental and calculated results for each particle size (d20, d60, d80,
and d90) when static equilibrium had been reached after 40 minutes of flow, shows that except
for the upstream reaches, the calculations generally agree with the experimental results. The
particularly good reproducibility (i.e., predictive accuracy) for d80 and d90 are likely a result of
the static equilibrium achieved through the entire waterway by that point. Thus, the authors’
technique is generally capable of explaining water surface profile, quantitative bed variation,and bed material size distribution in stony-bed rivers with a wide range of bed material sizes.

(a) After 10 minutes from the calculation start(a) After 10 minutes from the calculation start(b) After 30 minutes from the calculation start
Figure 9(b) After 30 minutes from the calculation startCalculated result of Sedimenttransport rate(c) After 40 minutes from the calculation start
Figure 8 Comparison between thecalculated result and the measured result(water surface and river bed elevation)

Figure 10 Comparison of the calculated
result and the measured result(D90, D80, D60, D20)

3.3 Differences in Fundamental Approaches between the Authors’ Technique andConventional Bed Variation Analysis Techniques
Applying conventional bed variation analysis to field experiments produces resultsthat differ greatly from the observed data. This section compares the fundamental approachesof conventional analysis and the authors’ technique and discusses problems with the former.
In general, conventional analysis uses sediment discharge equations, e.g., those of Ashida andMichiue (1972), to determine the transport of each particle size. Ashida and Michiue derive
their bed load equation for sediment mixtures using critical tractive force

Figure 11

Comparison of the tractive force and the critical tractive force

With sediment mixtures and relative share of each particle size in the bed materials inan equilibrium sediment discharge equation determined for uniform and fine sediment.
However, the sort of stony-bed rivers shown in Figure 1 contain large, greatly exposed,
immovable rocks, which prevent the sort of continuous bed material movement assumed by
Ashida and Michiue and instead result in the discontinuous movement of all bed material
sizes. The authors’ method accounts for such discontinuous movement and predicts thevolume of sediment discharge from bed material pick up volume and deposition volume, as inequations (8) and (9). Furthermore, bed load equation with mixtures derived by Ashida andMichiue uses critical tractive force with mixtures—an important indicator—based on amodified theory of Egiazaroff (1965). The Egiazaroff equation is:τ *c k =0.1{log 10 (19 d k d avg )}2(19)

In which critical particle transport velocity is calculated with a balance relation equation
assuming sliding motion. However, the manner of particle separation in a stony-bed river is
not sliding but, because of the extensive irregularities in the bed surface as shown in Figure 1,
is more likely explained by rolling motion in which one particle rolls over another. This is
why the authors make pick up judgment of particles in the manner seen in equation (4). In
Figure 11, tractive force at the point X = 40m after bed stability had been achieved as
calculated in Figure 8 is compared with the critical tractive force for each particle size as
calculated with the modified Egiazaroff equation. The modified Egiazaroff equation predicts
that, in a state in which the authors’ analysis predicts bed stability having been achieved,
particles smaller than 10 cm would be subjected to a tractive force exceeding the criticaltractive force and would therefore be transported.
Figure 11 also shows standard tractive force according to Parker (1990) (below which
only extremely little sediment discharge is believed to occur). Parker’s theory assumes thatparticles smaller than 5 cm are transported and so does capture the phenomenon to somedegree. However, because of the issue of porosity (discussed below) in conventional bed
variation analysis techniques, even if the concept of standard tractive force is incorporated,
the results will likely show a similar tendency as the Egiazaroff calculation results in Figure8(c).
Hirano’s equations (Hirano, 1971) are widely used as one-dimensional continuity
equations for mixed sand-and-gravel beds (only the equation for bed lowering is shown).

∂Z B
1∂ (q B ⋅ B )
a ∂λ
=−
+
(20)
∂t
B(1 − λ0 ) ∂X
1 − λ0 ∂t
The second term on the right side of the equation does account for temporal change inporosity. Because of the difficulty of calculating temporal change in porosity, however,porosity is usually treated as constant and the term is disregarded. Yet the large stones, etc., ina stony-bed river result in considerable surface irregularity and large voids around the rocks.
Calculations that treat the voids as constant indicate large rocks, which in reality would not
reduce in height, as lowering as the surrounding gravel is transported, yielding results
suggesting significant bed variation. By calculating height separately for each bed particle
size, the authors’ method accounts for porosity arising from bed irregularity, which is an
important mechanism determining bed variation in stony-bed rivers.

4.

CONCLUSIONS AND ISSUES FOR FUTURE STUDY(1) The authors have devised a new one-dimensional bed variation analysis method for stonybed rivers and validated this method by using it to reproduce results obtained in field
experimentation in the Joganji River. This validation showed a high degree of reproducibilityof important determinants of riverbed variation—the temporal and spatial distribution ofwater level profile, bed variation, and bed particle size distribution—and the method
reproduced the bed variation mechanism of stony-bed rivers to a degree not possible withconventional analysis.
(2) Further investigation is needed to determine the validity and range of applicability of thismethod in determining, for stony-bed rivers, the shielding effect on small bed particles oflarge bed particles, an important mechanism in this method.

REFERENCES
Ashida,K. and Michiue,M. (1972), Study on hydraulic resistance and bed-load transport rate
in alluvial streams, Proceedngs, JSCE, Vol.206, pp.59-69.
Egiazaroff, I.V. (1965), Calculation of Nonuniform Sediment Concentration, Proceedings,
ASCE, Vol.91, HY4, pp.225-247.
Fukuoka,S., Yamasaki,N., Kuroda,Y., Iuchi,T. and Watanabe,A., (2006), On river-bed
variation mechanism in a steep channel and method to calculate inundation discharge bydyke breach, Advance in River Engineering, JSCE, Vol.12, pp.55-60.
Fukuoka,S., Watanabe,A., Shinohara,Y., Yamasita,S. and Saito,K. (2005), Formation
mechanism of gravel bed rivers and it’s characteristics of grain size distribution, Advancein River Engineering, JSCE, Vol.11, pp.291-296.
Hirano,M. (1971), River-bed degradation with armoring, Proceedings, JSCE, Vol.195, pp.5565.
Kuroda,Y., Fukuoka,S., Yamamoto,T., Yoshida,K. and Iuchi,T. (2005), Formation
mechanism of gravel bed rivers and it’s characteristics of grain size distribution, Advancein River Engineering, JSCE, Vol.11, pp.363-368.
Osada,K. and Fukuoka,S. (2008), Development of 1-D river bed variation analysis method
based on sediment transport mechanism in stony bed rivers, Annual Journal of HydraulicEngineering, JSCE, Vol.51, pp.625-630
Parker.G. (1990), Surface-based Bedload Transport Relation for Gravel Rivers, Journal of
Hydraulic Research, Vol.28, No.4, pp.417-436.
Thorne,C.R., Bathurst,J.C. and Hey,R.D. (1987), Sediment Transport in Gravel-Bed Rivers,
Wiley.
Uchida,T., Fukuoka,S., Fukushima,T. and Tanaka,M. (2001), Two-dimensional shallow water
analysis over large roughness elements and its application, Journal of Hydraulic, Costal
and Environmental Engineering , JSCE, Vol.691, Ⅱ-57, pp.93-103.