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

Numerical analysis of flood flows and bed variations at riverconfluences of the Go River
Y. Takemura & S. Fukuoka

Research and Development Initiative, Chuo University, Kasuga, Bunkyo Word, Tokyo, Japan

ABSTRACT: In confluences of rivers having a comparable drainage area, a strong mixing of a mainstream and tributaries flows which have a large momentum in a different direction induces complex
three-dimensional flow structures. In this study, the Bottom Velocity Computation (BVC) method consistof shallow water equations, depth-integral vorticity equations and water surface velocity equations was
extended by introducing the flux integral method considering a non-orthogonal direction of a grid line forproperly evaluating flood flows and bed variations in the confluences. The extended model was applied torecent large floods in the Go River confluences to investigate the flood flows and bed variations during
the flood events. From the calculation results, we clarified that separation zones formed at the confluences
have pronounced effects on the three-dimensional flow structure and bed variations for the downstream
channel.
1

INTRODUCTION

The confluences of the Go River, Basen River and
Saijyo River are located in the Chugoku district,Japan. Each river has a comparable drainage area,
and the Basen River joins the Go River with agreat curvature at immediately downstream fromthe junction of the Basen River and Saijyo Riveras shown in Figure 1. In such river confluences, astrong mixing of the mainstream and tributaries
flows which have a large momentum in a different
direction induces complex three-dimensional flow
structures. Flood flows and bed variations at theconfluences are great concern of river engineers.

Figure 1. Air photograph of the Go River confluences(Sep. 2006).

Recently, full 3D turbulence model comes to be
used for investigating the three-dimensional flow
structures at river confluences (Miyawaki et al.
2010). But because the flow structures and sediment transport at the confluences are also affectedby spatial and temporal large scale flows, morepractical numerical model capable of evaluating theconfluence flows and bed variations is required.
In the case where the depth of a confluence area
is shallow, it is expected that three-dimensionalflows induced by deflections of depth-average flows
dominates in the flow structures and bed variations.
Recently, Uchida & Fukuoka (2011) proposed
the Bottom Velocity Computation (BVC) method
consist of shallow water equations, depth-integral
vorticity equations and water surface velocity
equations. Because the BVC method is based on
the shallow water model, three-dimensional flows
and corresponding sediment transports can be calculated more reasonably than the full 3-D turbulence model. On the other hand, since transports
and deformations of the depth-integral vorticity
are restricted by the depth-average flow fields, thecomputation accuracy of the depth-average velocity directly affects the resulting calculated threedimensional flow structures and bed variations. As
reported in previous studies, when the convection
terms in the momentum equations are approximated using upwind difference scheme, a large
error may result in the flow where the grid line
and flow directions are not closely aligned (Rhie &Chow 1983, Leonard et al. 1994). This influence
would be more pronounced in the confluences.

1811

In this study, we developed the flux integral
method considering a non-orthogonal direction of a grid line based on the CIP-CSL scheme(Nakamura et al. 2001) for evaluating the convection terms more accurately also in the complex
flow fields of the confluences and it was combinedwith the BVC method. The developed model wasapplied to recent large floods in October 1998 andJune 1999 at the Go river confluences to investigate the flood flows and bed variations during theflood events.
2
2.1

Since the convection terms of the depth-integral
equations of motion are calculated based on the
CIP-CSL scheme (Nakamura et al. 2001), the gridcell and grid boundary average unit-width discharge are solved by following equations:
s

s
1
∂ϕ
+ ( Fluxe − Fluxw + Fluxn − Fluxs ) = D (1)
S
∂t
l
l
l
l
l
l ∂ϕ
l ∂ϕ
l
l
∂ϕ
+U
+V
= −ϕ ∂U + ∂V  + D
∂t
∂x
∂y
∂y 
 ∂x

Computational method
Framework of computational method∂zs τ bx h ∂τ xx ∂τ xy Fx
+
−
+ 
−
− gh
∂y  ρρρ  ∂x∂xD=
∂zs τ by h ∂τ yx ∂τ yy Fy− gh ∂y − ρ + ρ  ∂x + ∂y  − ρ(2)Figure 2 shows framework of the computational
method of this study. The governing equations of
the Bottom Velocity Computation (BVC) method,i.e. shallow water equations, depth-integral vorticity equations and surface velocity equations(Uchida & Fukuoka 2011), are discretized on a­non-orthogonal grid to properly evaluate complex boundaries of river confluences. In the BVC
method, vertical distributions of horizontal velocity
components are approximated by a cubic function,and it allows calculating the momentum exchange
due to the vertical distribution of the velocity.

where ϕ (=Uh, Vh) = unit-width discharge components in x and y directions; U, V = depth-­averagevelocity components in x and y directions; h = water
depth; ρ = fluid density; g = the acceleration of­gravity; Fluxe,w,n,s = momentum flux from grid boundaries e, w, n and s shown in Figure 3. Here, overbars“–s” and “–l” indicate a grid cell and a grid ­boundary

Figure 2.

Figure 3. Arrangements of variables and interpolation
function for unit-width discharge distribution within agrid cell.

Framework of computation method.

1812

(3)

average value, respectively. Bed shear stresses τbx,y andvegetation resistances Fx,y are given as:τ bx = ρcbub ub2 + vb2 , τ by = ρcbvb ub2 + vb2
Fx =

ρ gha
ρ gha
U U 2 + V 2 , Fy =V U 2 +V 2
K
K

(4)
(5)

where ub, vb = bottom velocity components in x
and y directions; K = vegetation permeability coefficient; ha = height of a tree. cb is defined by the
relationship with the equivalent roughness ks:
cb = Ar +1  zb lnκ  ks (6)where Ar = 8.5; κ = 0.4; zb = δzb + ks; δzb = h/(e3 - 1).
In this study, ks is given as follows:
ks = (n/0.0417)

6(7)where n = the Manning’s roughness coefficient.
τxx, τxy, τyy, τyx are the horizontal shear stresses
consist of the momentum exchange due to theturbulence and vertical distribution of horizontal
velocity and they are written as follows:∂u+ u ′u ′∂x ∂v ∂u τ xy , τ yx = ν t ++ u ′v ′ ∂x ∂y ∂v
τ yy = 2ν t
+ v ′v ′∂yτ xx = 2ν t

(8)where vt = depth-average turbulent viscosity; u′,v′ =deviation velocity components from depth-­averagevelocity in x and y directions. vt is calculated fromequivalent condition of the turbulent kinetic
energy at each point. Here, overbar “–” indicates a
depth-average value.
Equations 1, 2 are solved in two steps. The firststep is to solve provisional values ϕ s* and ϕ l * fromthe equations without considering a right-handside of the equations (details are discussed in 2.2).
The second step is to correct the provisional values
by considering a right-hand side of the equations.
Through these steps, the values of ϕ s and ϕ l ofnext time step are obtained.
2.2

Computation method for convection terms

To calculate the convection terms of thedepth-­integral equations of motion accurately, the CIP_CSL method is introduced. The­two-dimensional CIP_CSL ­methods have been

proposed by researchers (Takizawa & Yabe 2000,
Nakamura et al. 2001). But because these methods
were proposed on an orthogonal grid, applications to the engineering problems were limited. In
this study, a calculation method of the convection
terms on a non-orthogonal grid will be developed
based on the CIP_CSL method.
Figure 3 shows the arrangements of the variables
used for the calculation of the convective terms. In
this model, distribution of ϕ within each grid cell isapproximated by a two-dimensional quadratic polynomial shown in Figure 3. 9 ­constants included in
the quadratic polynomial c0–8 are determined froma grid point values, a grid boundary average values
and a grid cell average value of ϕ at each grid cell(see Fig. 3). Here, the grid point values are interpolated from the grid boundary average values forthe advantage of the computation efficiency in thislstudy. A calculation method for ϕ and momentum flux at the gird boundary w is shown below
as an example. First, upstream points nw′ and sw′
described in Figure 3 are calculated from contravariant components of a depth-average velocity ϕ l atllgrid boundary w U cw , Vcw and the time step ∆t in
the general curvilinear coordinate. The provisional
value ϕ l* is obtained as the integral of the interpolation function of ϕ in grid cell W from nw′ to sw′.
The momentum flux from the grid boundary w arealso calculated by adding the integrals of the interpolation function in grid cell W throughout the areasurrounded by nw, nw′, sw, p and of the interpolation function in grid cell SW throughout the areal*surrounded by sw′, sw, p. ϕ and the momentum
flux at each grid boundary is calculated in the samemanner. And the provisional value ϕ S* of each gridcell is calculated from the budget of the momentum
flux from the grid boundaries e, w, n and s.
3

Application to the GO riverconfluences

Figure 4 shows the plane form and observation stations in the Go River confluences. The bed topography at the confluence consisted large bars, andthe excessive grows of vegetations and sedimentation have been essential problems. In this area,large floods occurred at October 1998 and June
1999 and flood data were obtained. We apply themodel constructed in preceding chapter to thesefloods for investigating applicability of the modeland characteristics of flood flows and sedimenttransports at the confluences.
3.1

Computation conditions

Observed water level hydrographs at Awaya(141.6 km point of the Go River), Minamihatashiki

1813

Figure 5. Grain size distributions in the Go River
confluences.

Figure 4. Plane form and observation stations in the
Go River confluences.
Table 1. The Manning’s roughness coefficient and vegetation permeability confident of calculated area.
Roughness coefficient(m-1/3 ⋅ s)
Go River (main channel)
Go River (flood plain)
Basen River (main channel)
Basen River (flood plain)
Saijyo River (main channel)
Saijyo River (flood plain)

0.032
0.042
0.032
0.045
0.032
0.042

*Vegetation permeability coefficients are given in therange from 30 to 50 at the vegetated areas.

(5.0 km point of the Basen River) and Miyoshi(1.0 km point of the Saijyo River) observation stations are given as upstream boundary conditionsof each river instead of discharge hydrographs.
Downstream boundary condition is given by aobserved water level hydrograph at Ozekiyama
observation station (139.0 km point of the GoRiver). The Manning’s roughness coefficient andvegetation permeability coefficient of calculatedarea are determined by reference to the previous
study (Fukuoka et al. 2007) as shown in Table 1.
Figure 5 is grain size distributions of the eachriver. Significant difference cannot be foundamong them. We give the grain size distributions
indicated by plots to each river uniformly throughout river reach for the calculation. In the calculationof sediment transport, bed load with 8 grain sizegroups ranging between 0.3 and 1.25 mm was considered. The bed load transport rate of each grain
size is calculated by Ashida & Michiue ­formula(1972). The change in the river bed ­elevation and

Figure 6. Comparison between observed and calculated discharge hydrographs of Oct. 1998 and Jun. 1999
floods.

fraction of the materials in surface layer are calculated by using Hirano equations (1971).
3.2

Computational results

Figures 6a, b are comparison between the observedand calculated discharge hydrographs of the1814

1998 and 1999 floods at the observation stations.
Although the calculated discharge hydrographs tend
to be large compared with the observed discharge
in the deciding period at Ozekiyama observation
station (139 km point of the Go River), they show
good agreement in the both floods. F­ igure 7 shows
longitudinal distributions of the calculated maximum water levels and flood marks of the 1998 and1999 floods. Although the calculated maximum
water levels almost reflect the height of the flood
marks, the difference between the calculated maximum water levels and flood marks becomes large
in the upstream reach of 2.0 km point of the BasenRiver. The main reason is considered below. The
Manning’s equation is used for resistance to flowin relatively low vegetation and the roughness coefficient is given uniformly for flood plains in eachriver in this study. As shown in Figure 5, the BasenRiver has relatively large flood plains, and distribution of the vegetation covering flood plains should
not be neglected for the resistant to flow.
3.3 Consideration for characteristics of flood
flows and bed variations in the Go River
confluences
Characteristics of flood flows and bed variations in the Go River confluences are investigatedbased on the observed data and calculation results.
­Figures 8a, b show calculated 10 minutes average

Figure 7. Comparison between the longitudinal distributions of the calculated maximum water levels and flood
marks in the 1998 and 1999 floods.

Figure 8. Comparison of distributions of 10 minutes average and instantaneous depth-integral vorticity in the direction of depth-average velocity at peak of the 1998 flood.

1815

and ­instantaneous distribution of depth-integral
vorticity in the direction of the depth-average velocity in the peak of the 1998 flood when the bed variations were occurred most. Red color corresponds
to the positive vorticity (counter-clockwise rotation)and blue color indicates the negative vorticity (clockwise rotation). The instantaneous vorticity distribution is quite different from the 10 minutes average
vorticity distribution, especially at the downstreamof the junction of the Basen River and Saijyo Riverand the flow convergence area around 139.2 km
point of the Go River, because of the unsteady vortex motion at the separation zones surrounded byblack broken line in Figure 8. In addition, the clockwise depth-integral vorticity (indicated in blue color)
formed at the downstream reach of the junction of
the Basen River and Saijyo River is transported
to the flow convergence area around 139.2 km pointof the Go River at the peak of the 1998 flood.
Figures 9a, b show the observed and calculated
bed variation contour before and after the 1998 and1999 floods. As shown in Figure 9a, deposition andbed scouring were occurred mainly at the junctionof the Basen River and Saijyo River and around theflow convergence area around 139.2 km point of the
Go River. The calculated bed variations shown in
Figure 9b can reproduce the bed scouring along theflow convergence area, but the scouring area is shortin a longitudinal direction and wide in a transverse
direction compared with the observed data shownby Figure 9a. In addition, the calculated bed variations cannot fully reproduce the depositions at thejunction of the Basen River and Saijyo River.
3.4

Consideration for effects of the momentum
flux of a non-orthogonal directionof a grid line

The key feature of the proposed model is considering the momentum flux of a ­non-orthogonal

direction of a grid line in the calculation of theconvection terms. Therefore, the calculation
without considering the momentum flux of anon-orthogonal direction of a grid line (Run-2) is
conducted under the same conditions as in the calculation of the preceding section (Run-1) to examine the applicability of the proposed model.
Figure 10 shows the calculated bed variation
contour before and after the 1998 and 1999 floods
by Run-2. It can be seen that the calculated bed
variations of Run-1 (Fig. 9b) and Run-2 (Fig. 10)­indicate different features at the downstream of theflow convergence area around 139.2 km point of the
Go River (139.2–139.0 km reach), and the calculation result of Run-1 comes closer to the observed
results indicated by Figure 9a than the result ofRun-2.
Next, trajectories of the water particles at vicinity of the water surface and river bed in the peak of
the 1998 flood are compared between Run-1 andRun-2 in Figure 11. The trajectories are estimatedfrom 10 minutes average calculated surface velocities
and bottom velocities. The trajectories at the vicinity of water surface assume a case observers toss
floating rods into the channel from the ­Kotobuki
and Tomoe Bridge shown in Figure 11. The trajectories at the vicinity of river bed are described byfocusing on the motion of the particles around theflow convergence area around 139.2 km. The discrepancy of the trajectories at vicinity of the water
surface between Run-1 and Run-2 becomes large
at the junction of the Go River and Basen Riverand downstream of the flow convergence area (seeFig. 11a). And the water particles at the vicinity ofriver bed toward the left bank in the downstreamof the flow convergence area in Run-1 comparedwith Run-2. It would be considered as the maincause of the difference of the calculated bed variations between Figure 9b and Figure 10 in thedownstream of the flow convergence area.

Figure 9. Comparison between the observed and calculated bed variation contour before and after the 1998 and1999 floods.

1816this, observations for the water surface velocities
and its trajectories are planned to be conductedby using air photographs and GPS floats during
floods. Based on these observed data, we advancethe examination for the applicability of the model.
4

Figure 10. Calculated bed variation contour before andafter the 1998 and 1999 floods (without considering a
non-orthogonal direction of a grid line in the calculationof convection terms).

Conclusions

This study developed the numerical model capable ofcalculating the flood flows and bed variations in theconfluences of rivers having a comparable drainage
area by introducing the flux integral method considering a non-orthogonal direction of a grid line on
the Bottom Velocity Computation (BVC) method.
The developed model was applied to the Oct.
1998 and Jun. 1999 floods of the Go river confluences, and its applicability and characteristics offlood flows and bed variations in the confluenceswere investigated together with the observed data.
From the investigation, we clarified that the separation zones formed at the junction of rivers andflow convergence area affects the three-dimensional flow structures of their downstream.
REFERENCES

Figure 11. Calculated trajectories of water particles atvicinity of water surface and river bed at the peak of the1998 flood in Run-1 and Run-2.

Although the applicability of the proposedmodel was investigated in this study, it is still inadequate because of the lack of observed data. The
proposed model can calculate trajectories of waterparticles at the vicinity of water surface which
can easily observed. In the Go River confluences,
water level gages have been set to obtain water­surface profiles in the confluences. In addition to

Ashida, K. & Michiue, M. 1972. Study on hydraulicresistance and bed-load transport rate in alluvial
streams, Proceedngs, JSCE, Vol. 206, pp. 59–69.
Fukuoka, S., Sato, H., Fujisawa, H. & Ohnuma, F. 2007.
Assessment of permeability coefficient of vegetation
and roughness coefficient under various floods andchannel conditions, Annual Journal of HydraulicEngineering, JSCE, Vol. 51, pp. 607–612.7.
Hirano, M. 1971. River-bed degradation with armoring,
Proceedings, JSCE, Vol. 195, pp. 55–65.
Leonard B.P., MacVean M.K. & Lock A.P. 1994. The
flux-integral method for multidimensional convection
and diffusion, NASA Tech. Memo., pp. 333–342.
Miyawaki, S., Constantinescu, G., Rhoads, B. &­Sukhodolov, A. 2010. Changes in three-dimensional
flow structure at a river confluence with changes inmomentum ratio, Proc. River Flow, Braunschweig,
Germany.
Nakamura, T., Tanaka, R., Yabe, T. & Takizawa, K. 2001.
Exactly Conservative Semi-Lagrangian Scheme forMulti-dimensional Hyperbolic Equations with Directional Splitting Technique, Journal of Computational
Physics, Vol. 174, pp. 171–207, 2001.
Rhie, C.M. & Chow, W.L. 1983. Numerical study of the
turbulent flow past an airfoil with trailing edge separation, AIAA Journal, Vol. 21, No.11, pp. 1525–1532.
Takizawa, K. & Yabe, T., 200. Development of Multi
Dimensional Conservative Semi-Lagrangian Scheme,
Symposium of numerical fluid dynamics, JSFM,
Vol. 14, pp. 1–6.
Uchida, T. & Fukuoka, S. 2011. A bottom velocity computation method for estimating bed variation in achannel with submerged groins, Journal of JSCE, Ser.
B1 (Hydraulic Engineering), Vol. 67, No.1, pp16–29,
in Japanese.

1817