Advances in River Sediment Research – Fukuoka et al. (eds)© 2013 Taylor & Francis Group, London, ISBN 978-1-138-00062-9

Study of flood flow and gravel riverbed variation analysis
in the Satsunai river
K. Osada

Anan National College of Technology, Tokushima, Japan

S. Fukuoka

Research and Development Initiative, Chuo University, Tokyo, Japan

H. Ohgushi

Obihiro River Office, Hokkaido Development Bureau, Hokkaido, Japan

ABSTRACT: In the Satsunai river with gravel-bed, stones and gravels have moved violently during
floods and severe bed scouring and bank erosion occurred. To clarify flow and sediment transport during
a flood is important to design a stable longitudinal and cross-sectional form. The authors have developed
a new method of two-dimensional riverbed variation analysis focused on bed surface unevenness and themechanism of sediment transport in gravel-bed rivers. We applied the new model and the conventionalmodel to riverbed variation of 2011 flood whose data of water surface profiles were observed in detail.
We showed that the new model is a useful method for flood flows and riverbed variation in the Satsunai
river with gravel-bed.
1

INTRODUCTION

In the Satsunai river with steep slope of theTokachi river system, stones and gravels are
moved violently during floods. Figure 1 showschannel changes over-time from 1971 average bed­elevation. Figure 2 shows cross-sectional form in
15 km and 19.2 km section. Channel form around1971 was wide and shallow. Many places of themain channel abutted on a river bank, because thealignment of the main channel meandered greatly.

A series of dikes have been constructed from 1978
for the bank protection of the main channel. But,
the main channel cross sectional form has transformed in narrow and deep shape. When a largescale flood happened, it is likely that severe erosionof flood plain and levee would occur. To clarify
flow and bed variation during a flood is important
to design a stable longitudinal and cross-sectional
form in the Satsunai river.
We develop a new model of two-dimensional riverbed variation focused on bed surface unevenness

Figure 1. Channel change from 1971s’ average bed­elevation (Dotted line: bank revetment, Black line: dike).

Figure 2.

589

Transition of riverbed cross-sectional form.

in gravel-bed rivers. In this model, to take accountfor the bed surface unevenness, we calculate heightsof each sediment size group on the bed surface.
Moreover, sediment discharge in gravel-bed rivers
is calculated by pick-up rate from the bed, deposit
rate to the bed and particle moving velocity. In
this study, we applied the authors’ model to 2011
flood for the check of validation, and compared
it with the conventional model which consists of
Hiranos’ sediment continuity equation (Hirano
1971) and Asida and Michiue sediment load equation ­(Ashida & Michiue 1972).
2

2.1

Two-Dimensional Modelof Riverbed Variation Analysis
in GRAVEL-Bed Rivers

Figure 4. Plan view and initial cross sectional form of
the field experimental channel.

Field investigation of bed surface unevenness

The authors (Osada & Fukuoka 2008, 2010)
expressed a height of each particle size by calculated average height of each particle size ZPk as
shown in Figure 3(a). However, we see to havevarious heights in the same size group such as
particles A and B in Figure 3(a). The differencesof these heights make the bed surface unevenness
and affect the amount of pick up and deposition.
Thus, the field investigation was necessary to elucidate the height distributions of each particle size
fzk based on the average heights of each particle
size ZPk as shown in Figure 3(b).
The field investigation was conducted in the fieldexperimental channel in the Jyoganji river (Maeshima
et al. 2011). Figure 4 shows a plan view and initial
cross-sectional form of the field-­experimental channel. The channel is 190 m long, 8 m wide and 1:130
bed gradient. The experiment was carried out insteady flow conditions (Case 1∼Case 4). The investigation was conducted on the bed surface after
Case 4. No. 0 and No. 16 sections in straight reachwere chosen as the investigation sections to ­eliminate

Figure 5.

direct influences of the channel meander and bankslope. Numbered seals were affixed on stones andgravel as shown in Figure 5, and heights of eachparticle ZPk were measured by surveying. Height ofa surface ZS is also measured to determine a baselevel of the bed. Each particle size was estimated by
the image analysis.
Figure 6 shows observed height distributions of
3 categories in particle size (d < 40 mm, 40 mm < d <75 mm, 75 mm < d). These profiles are close to thenormal distribution curve, and so equation (1) wasused for estimating each particle size distribution.
f ( z ′) =Figure 3. Relation between average height and height
distribution of each particle size.

Bed situation of investigation site.

 z′2 Z − ZPk1
exp  − 2  , z ′ = Pk
 (1)ZPk − Zs
2π σ k 2σ k 

Figure 7 shows observed height distributions
and normal distribution curves of equation (1).
­Figure 8 shows the relation between particle sizeand standard deviation of normal distribution
curves shown in Figure 5. Where dm = mean particle diameter. This relation is approximated byequation (2).

590

bed rivers. Figure 10 shows a procedure of the new
two-dimensional riverbed variation analysis. The
model is composed of unsteady two-dimensional
flood flow analysis and two-dimensional riverbed
variation analysis.
2.2.1 Unsteady 2-D flood flow analysis
The continuity and momentum equations of
unsteady two-dimensional flood flow analysis are
as follows:∂  h  ∂  q ξ  ∂  qη 
+
=0
 +
∂t  J  ∂ξ  J  ∂η  J 

Figure 6. Observed height distributions of each particle
size.

(3)

∂  q ξ  ∂  u ξ q ξ  ∂  uη q ξ ++∂t  J  ∂ξ  J  ∂η  J ∂ξ y q  ∂ξ∂ξ  q y  ∂ξ y+ uη− x  u ξ x + uη x  −  u ξ∂η ∂ξ∂η  J J ∂ξ22 ξx + ξ y ∂H ξxηx + ξ yη y ∂H 
= − gh 
+

∂ξ
J
∂η 
 J
Fξ
(4)
− D 90 + DξρJ
Figure 7.

∂  q η  ∂  u ξ qη  ∂  u η qη ++∂t  J  ∂ξ  J  ∂η  J ∂η y q  ∂η∂η  q y  ∂η y+ uη− x  u ξ x + uη x  −  u ξ∂η ∂ξ∂η  J J ∂ξ22 ξxηx + ξ yη y ∂H ηx + η y ∂H 
= − gh 
+

∂ξ
J
J
∂η 

Fη
(5)
− D 90 + DηρJApproximated normal distribution curves.

Figure 8. Relation between standard deviation and­particle size.


 d 

σ k = min 0.435 + exp  −6 k  , 0.6  
d

m



(2)

Maximum value of the standard ­deviation isassumed to be 0.6 because the value of ­standard
deviation do not become larger than 0.6.
­Equation (2) is used for estimating the bed surface
unevenness in the new model.
2.2 2-D riverbed variation analysis method
Figure 9 illustrates sediment transport mechanism
resulting from bed surface unevenness in gravel-where h = flow depth; qξ , qη = contravariant­components of discharge fluxes; uξ , uη = contravariant components of velocity vectors; H = waterlevel (h + ZB ); ZB = Average of riverbed elevation
J = Jacobian; ξx ,ξ y ,ηx ,η y = metrics; Dξ , Dη = theterms of depth-averaged Reynolds stress.
Form resistance of large stones is dominant to
bed resistance. The equations for evaluating bed
resistance are developed by using form resistanceof diameter d90 on the bed surface as follows:FDξ 90 = ξx Fx + ξ y Fy , FDη90 = ηx Fx + η y Fy 

(6) Fx 
 ux 
ε D 90
1
2
 F  = N D 90 2 ρCD AD 90uD 90 2
   (7)2 yux + uy  uy 591

Figure 9. Relationship between bed surface, flow and sediment transport.

Figure 10. Procedure of the two-dimensional riverbed variation analysis.

where Fx, Fy = x- and y-components of form
resistance of d90; N D90 = number of d90 on the bed
2
surface defined as N D90 = 0.2 α 2 d90; ε D90 = shielding coefficient of d90 (=0.48); ρ = water density;
CD = drag coefficient (=1.0); AD90 = projectedarea of d90; uD90 = flow velocity acting on d90
which is determined by the logarithmic velocity
distribution.

the particles (=π 4 , π 6 ). Each particle size fraction on the bed surface Pk is calculated by using
pick-up rate VPk and deposit rate VDk .
∂Pi , j ,k∂t

2.2.2 Average height and fraction of each sizeon the bed surface
Average height of each particle size ZPk is estimatedby using pick-up rate VPk , deposit rate VDk and eachparticle size fraction on the bed surface Pk .
∂ ZP i , j , k
∂t

(

)

α VP i , j ,k − VD i , j ,k=− 2
Pi , j ,kα3(8)where subscript i, j indicates the computational
grid number; α 2 ,α 3 = 2-D and 3-D shape factors of()

α 2 VP i , j ,k − VD i , j ,kdk
α3nk
α V i′, j ,k
+ P0i , j ,k ∑ 2
dk
k =1 α 3

=−

(

)

(9)

− VD i , j ,k, VPi , j ,k > VDi , j ,kVVi ′, j ,k =  P i , j ,k
, VPi , j ,k ≤ VDi , j ,k0where P0 i , j ,k = each particle size fraction in the­subsurface layer. Average of riverbed elevation ZB
is defined as equation (10).
nk

(

)

ZB i , j = ∑ Pi , j ,k ⋅ ZP i , j ,k −

592

k =1

dm i, j2(10)The first term on the right side represents theaverage particle height in the computational grid,and the second term is radius of the mean particle
diameter.
2.2.3 Estimation method of sediment transport
rate
Sediment transport rate of each particle size isevaluated by the product of the volume of sediment transport Vmk and particle velocity u pk .
qBξ i , j ,k = Vm i 0, j ,k u pξ i , j ,kqBη i , j ,k = Vm i , j 0,k u pη i , j ,k 

∂t+∂qBξ i , j ,k∂ξ+∂qBη i , j ,k∂ηConcept of the saltation analysis.

Figure 12.

Pick-up critical elevation.

Figure 13.

Rate of pick-up.

(11)

Volume of sediment transport per unit of area
Vmk is calculated by a balance of sediment transport rate and difference between pick-up rate anddeposit rate.
∂Vm i , j ,kFigure 11.

= VP i , j ,k − VD i , j ,k



(12)

Particle velocities of each particle size uPξ k , uPη k
are estimated by saltation analysis using the particle equation of motion as shown in Figure 11. Bed
surface is formed by the mean particle diameterconsidering the longitudinal and cross-sectional
bed gradient. Coefficient of restitution (=0.65)
is used for calculating the collision with the bed
materials.
2.2.4 Calculation method of pick-up rate
The particles located at high position such as
A particle in Figure 9 have a high probability ofpick-up from the bed. On the other hand, particles such as B and C particles cannot be pickedup from the bed by the effect of larger stones.
Pick-up equation was developed by considering
these mechanisms.
First, pick-up critical elevation of each particle
size is estimated by balance-of-moment equation(equation (13)). Pick-up critical elevation is definedas the elevation of β = 1 as shown in Figure 12, it iscalculated by shifting particle height.

α zkαρCDε aα 2d k2u 2fk + xk ρCLε aα 2d k2u 2fk
22
β i , j ,k =α xk ⋅ ( ρs − ρ ) gα 3d k3(β k ≥ 1.0 :Pick up, β k < 1.0 :No pick up)dd
α xk = k sinθ k , α zk = k cos θ k
2
2

(13)

where ε a = projected area coefficient (=0.7);
CL = lift coefficient (=0.85); ρs = particle density; u fk = the velocity acting on the ­particle;
θ k = angle of pick up from the bed. Next,
pick-up rate is determined by equation (14),when the sediment particle is judged to move by
equation (13).
VP i , j ,k = ε P′ i , j ,k ε ws i , j ,kN P i , j , kα 3d k3
TP i , j ,k



(14)

where ε Pk
′ = rate of pick up calculated by height­distribution; ε wsk = rate of pick up controlled by theshielding effect of large materials; N Pk = numberof particles of each particle size on the bed surface(= Pk α 2 d k2 ); TPk = time required for pick up fromthe bed.
Rate of pick-up of each particle size ε Pk
′ isdefined by a ratio above the pick-up critical elevation as shown in Figure 13, and it is calcurated byequation (15).

ε P′ i , j ,k = ε r ∫

∞

ZBk
′

593


1
z2 
exp  − 2  dz 
2π σ i , j ,k 2σ i , j ,k 

(15)where ε r: coefficient in consideration of a delay of
the time of the pick-up by the vibrancy.
ε wsk is estimated by equation (16).
 nk ′  Pw k ,k ′
Awk ′ ε wsk = 1 − max  ∑ 2 k ′=1  α 2d k ′nk ′  Ps k ,k ′+∑
Ask ′  , 1.0 
2
d
α


k ′=1
2 k′


3

3.1 Analysis condition(16)

where subscript k ′ indicates the particle size ofd k ′ ≥ d m and d k ′ ≥ d k ; Pw k ,k ′ , Ps k ,k ′ = respective
rates of particle k ′ located in a wake zone and at
the front of the particle k; and Awk ′ , Ask ′ = control
areas of pick up by shielding effects shown in­Figure 14 (α S = 1, αW = 1).
2.2.5 Calculation method of deposit rate
The amount of deposition on the bed per unit time
VD k is calculated by equation (17).
VD i , j ,k = PC i , j ,kVm i , j ,k 

(17)where, PC k is the rest ratio of each particle estimated by saltation analysis. As shown in Figure 15,
saltation analysis is conducted on the bed formed
so as to coincide approximately with particle size
distribution and height distributions, and moving
period ts is estimated by this saltation analysis. The
rest ratio PC k is calculated by equation (18) using
the moving period.
PC i , j ,k = ∫

1

1

0t
S i , j ,kt 
exp  −
 dt 
 tS i , j ,k 

Application to 2011 FLOOD(18)

We applied the new model and conventionalmodel (Egiazaroff 1965, Hirano 1971, Ashida &Michiue 1972) to the Satsunai river. The analysis
was performed from 14.0 km to 23.0 km as shownin Figure 16. The results of these models werecompared to the observed results of 2011 flood.
­Figure 16 shows observation stations of water levels. Water level gauges of simple pressure type were
installed at 6 points (No. 1~No. 6) in addition to
the existing 2 observation stations. The upstream
boundary condition was given observed waterlevel hydrographs at Daini-ohkawabashi (20.7 km
section) and No. 6, and the downstream boundary condition observed water level hydrographs at
Nantaibashi (15.0 km section).
Figure 17 shows air photographs taken before
and after 2011 flood. They show longitudinalchannel variations during the flood. To introduce
enough the initial riverbed form for the analysis isimportant for explaining these longitudinal channel variations. The reproduction of initial riverbed
form was made by referencing air photographs.
Figure 18 shows measured particle size distributions in the Satsunai river. These particle size
distributions were measured by removing stones
more than 0.15 m. Thus particle size distributionof black line shown in Figure 18 which was corrected was used for the analysis. Five particle sizes(250 mm, 120 mm, 50 mm, 10 mm, 2 mm) wereused in the calculation.
3.2 Result of analysis and consideration
3.2.1 On flood flow
Figure 19 and Figure 20 show comparison withobserved and calculated water level hydrographsand discharge hydrographs by using the conventional model. The Manning’s roughness coefficient

Figure 14.

Shielding effect by large stones.

Figure 15.
period.

Saltation analysis for estimated moving

Figure 16.

594

Observation stations of water levels.

Figure 17.

Channel variation during 2011 flood.

Figure 18. Initial grain size distributions for calculation.

n = 0.03 was given to reproduce water level and discharge for the conventional model. Figure 19 shows
the calculated results at No. 6 and Daini­ohkawabashi were generally reproduced observed
water level hydrographs. However, observed waterlevel hydrographs of the other points (No. 5, No. 3,
No. 1 and Nantaibashi) cannot be explained by theconventional model. The riverbed variations appear
to affect the water surface profiles. But, as shown
below, the result of riverbed variation by using theconventional model could not reproduce observed
riverbed variation, so that calculated water level
hydrographs were not reproduced observed results.
Calculated discharge hydrographs were estimatedexcessively in the falling period of the flood compared to the observed results.
Figure 21 and Figure 22 show the results of
the authors’ model. Figure 21 shows the calculated water level hydrographs at No. 6, No. 5 and­Nantaibashi were reproduced observed results
to a certain level. However, calculated water level
hydrographs at Daini-ohkawabasi and No. 3 were
estimated fairly low as compared to the observed
results. As shown below, the authors’ model cannot
explain enough bank erosion and channel variationof 2011 flood, and thus cannot reproduce waterlevel hydrographs. Figure 22 shows that calculated
discharge hydrographs agreed well observed results.
Figure 23 shows calculated results of depth andvelocity distribution at the time of the peak discharge in 2011 flood. However, the depth is deepand velocity is fast locally due to flow concentration to narrow and deep parts in the main channel.
3.2.2 Sediment transport and riverbed variation
Figure 24 shows the comparison between calculated and observed riverbed variation by 2011
flood. Figure 25 and Figure 26 show the calculated

Figure 19. Comparison between observed and calculated water level hydrograph (Conventional model).

595

Figure 20. Comparison between observed and calculated discharge hydrograph (Conventional model).

Figure 21. Comparison between observed and calculated water level hydrograph (Authors’ model).

Figure 22. Comparison between observed and calculated discharge hydrograph (Authors’ model).

Figure 23. Depth and velocity distribution at time of the peak discharge.

596

Figure 24.

Comparison between calculated and observed riverbed variation.

Figure 25.

Calculated result of D60.

Figure 26.

Sediment discharge of each grain size at time of the peak discharge.

results of D60 spatial distributions and sediment
discharge of each particle size.
First, the authors’ model result of riverbed
variation was under estimated longitudinally compared to observed riverbed variation as shown in
Figure 24. The data of sediment discharge by theconventional model indicate extensive transport ofgravel size materials as 10 mm, whereas larger size
materials were hardly transported. The calculated
sediment discharge and riverbed variation could

not estimate actual phenomena, so that calculated
results of water level hydrograph could not also
explain observed results. The conventional model
which cannot reproduce sediment transport oflarger size materials does not appear to explain riverbed variation in the gravel-bed rivers.
On the other hand, riverbed variation by authors’model was larger than the result of the conventional model as shown in Figure 24. ­However,
the authors’ model cannot reproduce enough597

bank ­erosion such as around Daini-­ohkawabashi(20.7 km) and bed scouring by the channel variation such as around No. 3 observation point,
because a mechanism of bank erosions are notconsidered. Figure 26 shows that large gravels
and cobbles were transported activity by using theauthors’ model, and its total transported amount
was larger than those of the conventional model.
The authors’ model cannot reproduce calculated
water level hydrograph in several points becausebank erosion and channel variation around theobservation point cannot explain enough. ­However,
the authors’ model is capable of explaining sedimenttransport rate and riverbed variation in the Satsunai
river compared to the conventional model.
4

CONCLUSION

We developed the new model of two-dimensional­riverbed variation focused on bed surface unevenness
and the mechanism of sediment transport in gravelbed rivers. We conducted the check of validation of
the authors’ new model and compared with the conventional model by 2011 Satsunai river flood.
The conventional model gives considerably lessresults than observed results of sediment transport rate and riverbed variation. The conventional model revealed to be insufficient to explain
sediment transport rate and riverbed variation in­gravel-bed rivers. On the other hand, the authors’model gives better sediment transport of large
gravel and cobble size groups and riverbed variation than the conventional model. This showedthat the authors’ model was a useful method for
sediment transport and riverbed variation in theSatsunai river with gravel-bed.
References
Ashida, K. and Michiue, M. 1972. Study on Hydraulic
Resistance and Bed-Load Transport Rate in Alluvial
Streams, Proceedings of the Japan Society of Civil
Engineers, No. 206, pp. 59–69.
Egiazaroff, I.V. 1965. Calculation of Nonuniform Sediment Concentration, Proc. ASCE, Vol. 91, HY4,
pp. 225–247.
Hirano, M. 1971. River-Bed Degradation with Armoring, Proceedings of the Japan Society of Civil Engineers, No. 195, pp. 55–65.
Maeshima, T., Iwasa, M., Osada, K. and Fukuoka, S.
2011. Study on River-Bed Variation and Its Grain
Size Distribution in Compound Stony-Bed River With
Meandering and Straight Channel, Annual Journal of
Hydraulic Engineering, JSCE, Vol. 55, pp. 769–774.
Osada, K. and Fukuoka, S. 2008. Development of Onedimensional Bed Variation Analysis Methods Focusedon the Mechanism of Sediment Transport in StonyBed Rivers, Advances in Hydro-Science and Engineering, Vol. 8, ICHE Nagoya, Japan.
Osada, K. and Fukuoka, S. 2010. Analysis of Two­Dimensional Riverbed Variation in Channels with
Stony Beds, 11th International Symposium on River
Sedimentation, South Africa, CD-ROM.

598