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

Three-dimensional gravel motions in numerical movable
bed channel with particles of various shapes and sizes
T. Fukuda, S. Fukuoka & T. Uchida

Research and Development Initiative, Chuo University, Bunkyo-ku, Tokyo, Japan

ABSTRACT: Sediment transport rate and bed variation in gravel-bed rivers are strongly affected bygravel sizes and shapes in streams. However, measurements of forces acting on sediment particles in active
are difficult, especially for rapid gravel-bed rivers. To investigate the dynamic mechanism of gravel moving
in streams, we have developed a numerical movable bed channel in which three-dimensional motions ofvarious shape and size gravel particles are calculated individually. The dynamic interaction between particles motions and hydrodynamic forces on movable bed surface in numerical experiment is investigated. It
is demonstrated that large gravel particles form static clusters and they have an important role to stabilize
gravel-bed surface against hydrodynamic forces of flows.
1

INTRODUCTION

The sediment transport process in gravel-bed ­rivers
during floods is explained that, at first small particles on the bed surface move, next large gravel
particles expose and resist mainly hydrodynamic
forces (Fukuoka et al. 2007, Osada and Fukuoka2012). Exposed large particles do not move ­easily
in large floods, rest for a long period, and rollintermittently. However, derivation of sediment
discharge formulae for sand assumed that all particles move to a depth of several particle diameters.
This sediment transport mechanism differs fromthat of gravel-bed rivers. Bed variation analysis
using sand discharge formulae can not evaluateadequately the sheltering effect of large particles
on small particles (Suzuki et al. 1994) and hydrodynamic forces acting on projecting large particles.
Therefore, it has been pointed out that an analysis
of bed variation in gravel-bed rivers using thesesediment discharge formulae is questioned (Osadaand Fukuoka 2012).
Many researchers have investigated the sediment transport mechanism in gravel-bed rivers inwhich large gravel particles move intermittently,
and reported that formations of large particles
clusters on the bed surfaces contribute to river
beds ­stabilities. Piedra (2012) made experiments
for movable bed of graded sediments and investigated the relationship between tractive forces andratio of large gravel cluster occupying area to the
bed surfaces. Hendrick (2010) made field observation and reported distribution of large particles
clusters and their changes due to flood events.
Characteristic structural arrangements of particles

in gravel-bed rivers, such as imbrication, indicatesthat stabilities of river beds are strongly affectedby the sizes and shapes of particles. However,experiments and field observations have difficulty
as ways to investigate the dynamic mechanism ofgravel transport and how particles are arranged
against hydrodynamic forces in gravel-bed rivers inwhich bed materials are widely distributed in sizes
and shapes.
In recent years, progress of computing performances have made possible to analyse the transportmechanism of individual particles motions by the
Lagrangian method (Goto et al. 2000). In most ofthese analyses, hydrodynamic forces acting on particles are calculated using drag coefficient therefore
evaluations of hydrodynamic forces have remained
inexact. Ushijima et al. (2008) developed numericalmethod for solid-liquid multiphase flow in whichhydrodynamic forces are calculated directly by calculating three-dimensional fluid motion around asolid object using finer computational cells thanthe object size. Harada et al. (2011) conductedlarge eddy simulation of solid-liquid multi-phaseflow for drifting particles on movable bed of threesizes spheres and investigated vertical sorting process in oscillatory flows. Fukuda et al. (2012) has
developed a numerical model evaluating threedimensional gravel particle motions in consideration of various shapes of particles by applying thenumerical method of solid-liquid multiphase flow
to solid-liquid interactions and combined spheres
model (Matsushima et al. 2009) to various shapeparticles (Fig. 1), and demonstrated that gravel
particle trajectories with frequent saltation andlateral migration in experiments could not be323

Figure 2.

Figure 1. Gravel particles model.

r­ eproduced in the simulation unless the irregularshapes of particles are considered. Using the developed model, a numerical movable bed channel withparticles of various shapes and sizes is constructed
and mechanism of gravel transport and dynamic
effect of gravels to the streams are investigatedusing the numerical channel.
2

NUMERICAL MODEL

Concept of the fluid analysis.

where ui: i-direction averaged velocity given weightby the mass within a fluid cell, P: total of pressure
and isotropic component of SGS stress, ρ: density,µ: dynamic viscosity, gi: gravitational acceleration,
vt: SGS turbulent viscosity, CS: Smagorinsky constant ( 0.173 ), ∆: calculation grid size. ­Physicalproperty φ (density ρ, dynamic viscosity µ) andaveraged velocity ui are calculated as volume­averaged values and mass-averaged values, respectively, as follows:

In the present numerical model, fluid motions
are simulated with the Eaulerian method andgravel motions are simulated with the Lagrangian
method. To take into account the effect of solid
phase on liquid phase, fluid motions are simulated by governing equations of multiphase flows(Fig. 2). Fluid forces acting on particles are evaluated directly by integrating the forces acting on
particles in the multiphase flow. Gravel particles
with various shapes are made by the combinationsof several small spheres (Fig. 1). The motions of
the particles are simulated as rigid bodies.

φ = αφs + (1 − α ) φ f , φ f = f φl + ( 1 − f ) φ g

2.1

2.2 Evaluation of solid rigidity in fluid calculation

Governing equation of fluid motion

Fluid motion is simulated by one-fluid model for
multiphase flow using Smagorinsky model as the
subgrid turbulence model:
∂ui
=0
∂xi

(1)

∂u
∂ui
∂
1 ∂P
+
+ u j i = gi −2 (ν + ν t ) Sij∂x j∂tρ ∂xi ∂x j

{

Sij =

1  ∂ui ∂u j +
2  ∂x j ∂xi 

ν =µ ρ
ν t = ( CS ∆ )

}

(2)
(3)
(4)

2

2Sij Sij(5){

ui = αρs usi + (1 − α ) ρ f u fi

}ρ

(6)
(7)

where f: fluid volume fraction in fluid calculation
cell, α: solid volume fraction in fluid calculation
cell, suffix l, s and g denote liquid phase, solid
phase and gas phase, respectively. uf is fluid velocity and us is solid velocity. Continuity equationand momentum equation are solved by the SMAC
scheme with staggered grid.

Solid rigidity is evaluated by simulating particle
motion as rigid body and setting average density
ρ and average velocity ui with eq. (6) and eq. (7).
Solid velocity usi in eq. (7) is calculated equation
below:
us = rG + ω × rf

(8)

where rG : translational velocity vector of particle
gravity center, ω: angular velocity vector of particle, rf: position vector from gravity center to theevaluation point of us. The solid volume fractionsin fluid calculation cell α in eq. (6) and eq. (7) arecalculated by splitting a fluid calculation cell into
sub-cells and counting the number of sub-cells
included in gravel particle (Fig. 3).

324

2.5 Governing equation of particle
motion in streams
Particle motions are simulated by momentum
equations and angular momentum equations as
rigid bodies:MrG = Mg + Ff + Fc

{ (

(10)

)

.
ω r = I r −1 R −1 N f + Nc − ω r × Irω r
Figure 3.

Evaluation of particles volume with sub-cell.

Figure 4. Cells image for calculation of rigid bodyproperties.

2.3

Evaluation of water surface variation

Water surface variation is evaluated with continuity equation of fluid volume fraction f based on the
VOF method (Hirt and Nichols, 1981).
∂f ∂fui=0+∂t ∂xi(9)

Gas motion is not calculated to reduce calculation time and fluid motion is calculated underpressure zero condition at the water surface asboundary condition.
2.4

Modeling of various shape particles

Various particles in shape and size are made compactly by the combination of several small spheres
to detect contact and collision points between particles as seen in Figure 1. It is difficult to calculate
shapes of lapped part of spheres geometrically,
therefore, evaluation of rigid body properties suchas mass, gravity center and tensor of momentum
inertia become major issues. We calculate rigid bodyproperties of particles by using sufficiently small
cells for calculation of rigid body properties andnumerical integration in consideration of volumes
and points of cells included in particles (Fig. 4). The
rigid body properties of particles are not changed
in motion, therefore they have only to be calculated
just once at the beginning of the simulation.

}

(11)

where bold face letters indicate vector tensor andmatrix, M: mass of particle, rG: gravity center
point, g: gravitational acceleration, F: forces
acting on particles surfaces, N: torques acting
on particles, ω: angular velocity, R: transformation matrix from the local coordinate system tothe global coordinate system, and I: tensor of
momentum inertia. Suffix f, c: fluid force andcontact force, respectively, suffix r:components
in local coordinate systems of each particle. Dot
notations over letters denote time derivatives. In
calculation of angular momentum equations, atfirst, next time step angular velocities of local
coordinate system ωr are solved with eq. (7), next,
angular velocities of global coordinate system ωare calculated with transformation matrices R.
Next time step transformation matrices R are set
in consideration of the rotation in time step ∆twith the angular ­velocities ω. Next time step small
sphere positions are set with next time step transformation matrices R and small sphere positions
in local coordinate system. In these coordinatetransformation between global coordinate system
and local coordinate system, quaternion is usedinstead of transformation matrix R (Ushijima
et al. 2008).
2.6 Evaluation of fluid force
Fluid forces acting on a particle is evaluated withintegrating forces acting on a particle in the multiphase flow (Ushijima et al. 2008):
F f ,i = ∫

 ∂P
∂
+ρ
2 ( ν + ν t ) Sij
−
Ωs
∂
x
∂
xj
i


{



}  d Ω (12)



∂
 ∂P
N f ,i = ∫ ε ijk rf , j −+ρ
2 (ν + ν t ) Skl  d Ω
Ωs
∂xl

 ∂xk(13)
where Ff,i: i component of fluid force, Nf,i: i component of torque caused by fluid force, rf,i: position
vector from the gravity center of the particle to thefluid calculation cells, Ωs: an area occupied by aparticle and ε ijk Levi-Civita symbol.

325

{

}

Figure 6. Shape of the numerical movable bed channeland definition of coordinate axes.

Figure 5. Evaluation of contact forces.

2.7

Evaluation of contact force

Contact forces acting between particles are calculated by the contact detection of each small spherecomposing a particle (Fig. 5). Contact forces andtorques acting on gravity centers of particles arecalculated by the summation of the contact forces
calculated concerning each small sphere:
Fc = ∑ Fcp,n Nc = ∑ rcp,n × Fcp,n,(14)

where, Fc, Nc: contact force and torque acting on
the gravity center of a particle, Fcp,n : contact force
acting on each sphere, rcp,n : position vector fromthe particle gravity center to the contact point.
Contact forces between particles are calculated by the Distinct Elements Method (­Cundall
and Strack, 1979) and the spring constants k
and coefficients of dashpot c are calculated by­equations (15)–(17) (Gotoh, 2004):
1

2
 4  r r   3E 
kn =   1 2  
e

n

2
 9  r1 + r2   1 − pos 


s0 =

ks
1
=
kn 2( 1 + pos )cn = 2 h

m1m2
kn , cs = cn s0
m1 + m2

(15)
(16)
(17)

where suffix n and s: a component of direction
from center of spheres to the contact point andorthogonal two directions on the tangential plane,
respectively, r1, r2: radii of contacting two spheres
m1, m2: mass of contacting spheres.
3

NUMERICAL EXPERIMENTON MOVABLE GRAVEL BED
WITH various shape and size

3.1

The out line of numerical movable bed
channel and experimental condition

Figure 7.

Grain size distributions.

numerical channel with 15 m long 1 m wide and1:20 bed gradient was designed by considering flow
condition entraining gravel particles and length of
bed waves and computational load. Four different
shapes of particles are prepared which are composedof 8–10 small spheres as shown in Figure 1. Fivesizes of particles (40 mm, 50 mm, 70 mm, 90 mm,
120 mm) are prepared by varying small spheres
sizes composing particles. The diameters of thevarious shape particle are given as the diameter ofspheres with the same volume. Particles are packed
in to the channel by varying the supplying rate
to fit the particle size distribution in Figure 7. In
this paper we define large particles as larger than
90 mm and small particles as smaller than 70 mm.
The discharge of 0.5 m3/s is given at the upstream
end of the channel, and zero pressure condition is
given at the downstream end. The number of particles discharging out the end of the channel aresupplied simultaneously at the upstream point in
x = 1–2 m section. The numerical analysis of the
movable bed was started at the same time as thestart of discharging water, and we defined the timeof numerical experiment start as water discharge atthe downstream end becomes almost constant, andthe numerical experiment was conducted in 70 s.
Hereinafter experimental time is defined as the timefrom experiment start. The parameters used in thesimulation are shown in Table 1.
3.2 The result of the numerical experiment

Figure 6 shows the shape of the numerical channel with movable bed and definition of coordinate
axes used for numerical analysis. The shape of the

3.2.1 The outline of the variation of the flow
condition, bed condition and sedimenttransport discharge
Figure 8 shows the time variation of longitudinal
bed levels and water surface profiles averaged in

326

Table 1.

Parameters used in the simulation.

Δx, Δy, Δz: Fluid calculation cell size
Δt: Time Step for fluid calclation
ρw: Density of waterρs: Density of gravel
μw: Dynamic viscosity of water
μs: Dynamic viscosity of gravelΔt′: Time step for DEM calclation
E: elastic modulus
pos: Poisson’s ratio
h: Coefficient of dashpot
Number of gravel particles
Number of small spheres
composing gravel particles

0.01 (m)5.0 × 10-4 (s)1,000 (kg/m3)
2,650 (kg/m3)
8.9 × 10-4 (Pa ⋅ s)
8.9 × 10-4 (Pa ⋅ s)

Figure 8.

Longitudinal profiles of averaged water levels and bed levels.

Figure 9.

Space and time variation of bed surface.

lateral direction. The values are averaged for two
seconds. After 70 s from experiment start, average
depth is 0.23 m, average velocity is 2.17 m/s andFroude number is 1.45. The flow is supercriticalflow and water surfaces vary in the same phase as
the bed wave. Figure 9 shows numerical computation results of the bed surface. Bed level and distribution of particles diameter on the bed surface
are shown in Figure 10. In Figure 9 and Figure 10,
the bed condition before the experiment is also
shown for the reference. Regarding the bed profiles
shown in Figure 8, bed degradation occurs from

5.0 × 10–6 (s)5.0 × 1010 (Pa)
0.33 (–)
0.11 (–)
33 × 103 (–)
301 × 103 (–)

the ­downstream end, and it propagates toward
upstream and bed surface wave comes to be seen.
Relatively high positions are seen at x = 4 m, 7 m,
9 m and 13 m points in Figure 8. At those points
large particles form clusters and particles sorting
occurs on the bed surface. Figure 11 shows sediment discharge rate of each particle size measured
at the downstream end, and Figure 12 shows thenumber of flowing out particles from downstream
end per unit time. The sediment discharge rates
shown in Figure 11 vary in short time although
total transport rate of each particle is almost327

Figure 10. Bed level and particle diameters distributions on the bed surface.

Figure 11. Sediment discharge rates.

Figure 13. Snap shot view of the flow vectors near the
bed surface (t = 45 s right side of the channel).

Figure 12. Number of flowing out particles.

s­ table. Sediment discharge rate of large particles is larger while flowing out numbers of small
particles are larger in Figure 12, therefore it wasrecognized that small particles are also moving
well as large particles. In this experiment, aggradations occur at supplying section x = 1–2 m andx = 4–5 m, and degradations in the downstreamof x = 5 m point, and bed is not in equilibrium
condition. However sediment discharge rate isalmost stable and bed level is also stable in downstream of x = 5 m point, while particles overall
sizes are moving. Therefore it is ­considered that
this result of numerical movable bed experiment isworthy to investigate movement and arrangementof gravel particles against streams and mechanismof particle motions.

Figure 14. Vertical distribution of flow and particles
x-direction velocities.

3.2.2 Factor of increasing transport
rate of large particles
The factors of large amount of the sediment discharge of large particles, larger than 90 mm, areinvestigated. Figure 13 shows vectors of streams
near the bed surface, and Figure 14 shows the ­vertical
distributions of time averaged flow ­velocities in thex-direction at the center of the channel. The velocities are averaged over 0.15 m lateral range which isabout the largest particle size. Figure 14 also showstime averaged velocities of particles passing thesection. As seen in Figure 13 and F­ igure 14, fluid

328

velocities decrease considerably from upper position
to the lower position within the scale of larger particles scale and is about 2.0 m/s at the upper positionof large particles while at the lower position of thelarge particles where small particles are arranged,
fluid velocity decreases greatly and is smaller than
0.5 m/s. The large particles located relatively highposition received large velocities, therefore sedimenttransport rate of large particles become great.
3.2.3 The process of large particles
cluster formation
We investigate the factors of large particles cluster
formation on the bed surface shown in Figure 9
and Figure 10. Figure 15 shows trajectories of eachsize particles from t = 42.5 s to t = 47.5 s. In addition to the figures of the x-y planer trajectories,
bed level is also shown on the background. The trajectories of small particles, smaller than 70 mm, in
x-y plane take deep water passage beside projecting beds, while large particles, larger than 90 mm,
move over projecting beds formed at center of thechannel around x = 8.5 m point and right side andleft side of the channel around x = 13 m point. It
is because gravity centers of rolling large particles
are high position. Therefore large particles receive
upward force due to collision with projecting beds
at lower part than gravity center of particles, andsmall particles cannot move over projecting beds

Figure 15.

but come to take detours because gravity centerof small particles are low position on the bed surface. In this manner, large particles move over projecting beds more easily than small particles andreceived large fluid velocity, therefore large particles are not able to stop easily on flat bed. Largeparticles can stop at projecting beds as high aslarge particles size, therefore it is considered thatlarge particles tend to form clusters. Deeper water
passages are formed beside large particles clustersof the bed surface and small particles come tomove deeper water passages beside large particles
cluster. In this manner, particles on the bed surface are divided into large particles clusters andother particle groups. This results in the particle
sorting of the bed surface.
3.2.4 Hydrodynamic forces acting
on the bed surface
To investigate how particles resist against streams,
we convert the hydrodynamic forces acting on particles to the forces acting on unit x-y plane area,
evaluating hydrodynamic forces as the stress, andinvestigate the relationship between converted
stress, uneven bed surface shape, and sizes ofparticles on the bed surface. Figure 16 shows­longitudinal profile of bed level in upper figure,longitudinal profiles of stress acting on particles
in middle figure, and longitudinal profiles of mean

Trajectories of particles (t = 42.5 s–47.5 s).

Figure 16. Longitudinal profiles of water level, bed level, hydrodynamic forces acting on particles, mean diameter ofparticles on the bed surface (t = 45 s).

329

particle sizes on the bed surface in lower figure.
The stress acting on particles from streams areevaluated with eq. (18) by integrating x-directionhydrodynamic forces from channel bottom to the
top of the moving particles and dividing the areaof x-y integrating area to evaluate forces acting on
per unit bed surface area.

τ f , zx
=


 ∂P

∂



∫ ∫ ∫ α − ∂x1 + ρ ∂x j {2 (ν + νt ) S1 j } dxdydz
∫ ∫ dxdy

(18)

Integrating range in lateral direction is fromy = −0.15 m to y = 0.15 m and longitudinal integrating distance is 0.01 m as fluid calculation cells
sizes. Mean particles sizes in lower figure indicate
the averaged value of 0.15 m squared area on the
bed surface evaluated as area-averaged diameter ofparticles seen on the bed surface. From Figure 16,
it was demonstrated that resisting forces against

streams were acting at the projecting beds formedby large particles clusters.
3.2.5 The effect of particles shapes on startingand stopping of their motion
We investigate how particles behave when particles
start to move and stop. To investigate the behaviorof particles, we define the eigenvector corresponding to the largest value of principal moment ofinertia as short axis, and eigenvector corresponding
to the shortest value as long axis. A short axis anda long axis defined in this manner are orthogonal.
A long axis is the easiest axis to rotate among principal axes. Short axis is the hardest axis to rotate,
and indicates almost the same direction as normal
vector of the most flat plane of particles. Longaxes and short axes of the particles are shown in
Figure 17. Figure 18 shows time variations of averaged particles velocities in x-direction, and contact
forces and hydrodynamic forces acting on particles
nondimensionalized by x-direction gravity force,and angles between a long axis and x-axis andangles between short axis and z-axis, the time ofthe data was selected concerning the starting and

Figure 17. Conditions of long axis and short axis of particles on the bed surface.

Figure 18. Time variation of the particles velocities, hydrodynamic forces, contact forces and particles direction­concerning starting and stopping motion.

330

stopping of particles motion, and all the values areevaluated by averaging individual particles motion.
In this figure, red lines in upper figures show particles velocities and starting of the particles motion
is shown as a time detaching from zero, and stopping of the motion is a time corresponding to zero.
The time zero is selected at the time when a particle
begins to move and comes to a stop. In the figureof starting motion of large particles, it is seen that,at the starting motion of large particles, the angle
between long axis and x-axis are becoming larger.
According to the variation, it is revealed that static
large particles incline long axes toward downward,
and at the beginning of the motion, they change
the direction of long axes to roll downward easily. The angle between short axes and z-axis also
become larger, thus static large particles point flat
plane toward up direction, and at the starting, theychange the directions of flat plane. Accordingly
projected area is increasing, and therefore hydrodynamic forces shown in upper figure is increasing
together with the variation of direction of ­particles.
This kind of characteristics of the time variationof the motion are shown in the figure of small
particles, however, characteristics of the time variation is not seen so clearly. This is because that, atthe lower part of the bed surface where small particles rest on, stream direction changes largely due toarrangement of the static large particles and small
particles rest in a gap of large particles, therefore,
small particles starting and stopping motions areaffected by the irregular shape of gap. Concerning
the stopping of large particles, time variations oflong axes and short axes show the opposite variations, long axis inclines to the x-direction andshort axis to the z-direction then it was revealed
that particles end up the motion by inclining theirflat planes in the upward direction. Characteristic
variations are that the large particles motion are
stopped by receiving negative x-direction contact
forces and accordingly hydrodynamic forces acting
on the particles are increasing a little, on the other
hand, hydrodynamic forces acting on small particles are decreasing after the stop of small particles
motions. The reason of the differences is considered that small particles can stop at small velocity regions behind stationary large particles, whilelarge particles stand high positions and receive
large hydrodynamic forces and therefore cannot
stop easily but do by the collision and support ofsurrounding large particles.
4 cONclusion
This paper developed a numerical movable bed
channel composed of particles with various shapes
and sizes and predicted the three-dimensional­ articles motions and forces acting on streams
pfrom particles which are difficult to measure in
gravel-bed rivers. The results of the numericalexperiment draw following principal conclusions.
1. From the investigation of detailed particles
motion data, the motion mechanism of eachsize particles in streams and resistant mechanism of different sizes and shapes particles instreams were demonstrated.
2. Large particles start to move with inclining the flat
plane direction from upward direction to anotherdirection, accordingly hydrodynamic forces are
increasing, and inclining the long axis from downward to another direction to roll easily. Small particles indicate similar time variation although thetrend of the variation is not so clearly than large
particles, because small particles are affected byirregular shapes of gaps formed by surrounding
static large particles and streams around small
particles are disturbed and assume various direction due to surrounding large particles.
3. Rolling and saltating large particles move highposition of the bed surface and receive large flow
velocity, therefore large particles cannot stopeasily. Large particles stop by collision with surrounding large particles and tend to form large
particles clusters. Small particles move deeper
water passage beside clusters therefore particle
sorting occurs at the bed surface. The formationof large particles clusters has an important roleto stabilize gravel-bed surface resisting against
large hydrodynamic forces of flows.
References
Cundall, P.A. and Strack, O.D.L. 1979. A discrete numerical model for granular assemblies, Geotechnique,Vol. 29, No. 1, pp. 47–65.
Fukuoka, S. and Abe, T. 2007. Mechanism of low water
channel formation and the role of grain-size distribution in gravel-bed rivers, 10th International Symposium
on River Sedimentation, Moscow, Russia.
Fukuda, T., Fukuoka, S. and Uchida, T. 2012. Threedimensional numerical modeling for the motion ofstones with different sizes and shapes in streams, Proc.
tenth International Conference on Hydro-Science andEngineering, Orlando, USA.
Gotoh, H. 2004. Computational mechanics of sedimenttransport, Morikita Shuppan Co., Ltd, in Japanese.
Gotoh, H., Yeganeh, A. and Sakai, T. 2000. Coupling ofmultiphase-flow model and distinct element method
for simulation of sediment transport under high bottom shear, Journal of Hydraulic, Costal and Environmental Engineering, JSCE, Vol. 649/II-51, pp. 17–26,
in Japanese.
Harada, E., Tsuruta, N. and Gotoh, 2011. Large eddy
simulation for vertical sorting process in sheet-flow
regime, Journal of JSCE, B2 (Costal Engineering),
Vol. 67, No.2, pp. I_471–I475, in Japanese.

331Hendrick, R.R., Ely, L.L. and Papanicolaou, A.N. 2010.
The role of hydrologic processes and geomorphology
on the morphology and evolution of sediment clusters in gravel-bed rivers, Geomorphology, Vol. 114,
pp. 483–496.
Matsushima, T., Katagiri, J., Uesugi, K., Tsuchiyama, A.
and Nakano, T. 2009. 3D shape characterization andimage-based DEM simulation of the lunar soil simulant FJS-1 J. Aerospace Engineering, Vol. 22, No. 1,
pp. 15–23.
Osada, K. and Fukuoka, S. 2012. Two-dimensionalriverbed variation analysis method focused on themechanism of sediment transport and the bed surface unevenness in stony-bed rivers, Journal of JSCE,B1 (Hydraulic Engineering), Vol. 13, pp. 339–344, inJapanese.

Piedra, M.M., Haynes, H. and Hoey, T.B. 2012. The
spatial distribution of coarse surface grains and the­stability of gravel river beds, Sedimentlogy, Vol. 59,
pp. 1014–1029.
Suzuki, K., Kurihara, T. and Yamamoto, H. 1994. Studyon mixing effects on sediment discharge of sand andgravel with uniform grain size, Journal of Hydraulic,
Costal and Environmental Engineering, JSCE, Vol. 38,
pp. 627–632, in Japanese.
Ushijima, S., Fukutani, A. and Makino, O. 2008.
­Prediction method for movements with collisions ofarbitrary-shaped objects in 3D free-surface flows,
JSCE Journal B, Vol. 64 No. 2 pp. 128–138, inJapanese.

332