Numerical analysis of unsteady flow and bed variation using temporalchanges in water surface profiles during 1981 flood of the Ishikari river
mouth
Seiji Okamura1 and Shoji Fukuoka2
1

Department of Science and Engineering,
Chuo University, Japan
E-mail: seijiokamura5@gmail.com
2
Research and Development Initiative,
Chuo University, Japan
E-mail: sfuku@tamacc.chuo-u.ac.jp

Abstract
For the flood control and the river management, it is important to evaluate the degree of bed
variation during floods. But accurate measurements of temporal changes in bed elevations
during floods are difficult. In this study, we develop the quasi-3D unsteady flood flow and 2D
bed variation analysis method using observed temporal changes in water surface profiles. This
computational method utilizes the fact that the influences of channel shape, bed resistance andbed variation are reflected in temporal changes in water surface profiles. In the present
method, the concentration of suspended sediment is calculated by 3D advection-diffusion
equations. This paper shows that the present computation method can explain properly bed
variation mechanism during the August 1981 flood of the Ishikari River which a large bed
scouring at the river mouth occurred.
1.

Introduction

For the flood control and the river management, it is important to evaluate the degree of bed
variation during floods. But, it is difficult to elucidate temporal changes in bed elevations
during the period of floods because of the difficulty of measurement.
Fukuoka et al. (2004) estimated the bed resistances and the discharge hydrographs bycomputing 2D unsteady flood flow so as to coincide with observed temporal changes in water
surface profiles. This computational method is based on the idea that the influences ofchannel shape and resistance are reflected in temporal changes in water surface profiles. The
influence of bed variation during flood is also reflected in temporal changes in water surface
profiles where the large degree of bed variation occurs. For the purpose of river management,
Fukuoka (2011) proposed an observation-calculation system of flood flow and sediment load
focusing on temporal changes in water surface profiles to figure out flood flows and bed
variations during floods. In the system, temporal changes in water surface profiles anddischarge hydrographs during the period of floods and bed elevations before and after floods
are observed at the sections in question. Then, we can figure out flood phenomena bycomputing the unsteady flood flow and bed variation so as to explain these observed data.

1

The objective of this study is to demonstrate that the bed variation mechanism during theAugust 1981 flood of the Ishikari River which a large bed scouring occurred in river mouth
can be evaluated by the quasi-3D unsteady flow and 2D bed variation analysis using observed
temporal changes in water surface profiles. In the analysis, the concentration of suspended
sediment is calculated by 3D advection-diffusion equations in order to evaluate suspended
load affected by the secondary flow in meandering channel.
2.

The August 1981 Flood of the Ishikari River Mouth and the Previous Researches

2.1 Observation at the August 1981 Flood of the Ishikari River Mouth

Figure 1 shows observed discharge hydrograph of August 1981 flood of the Ishikari Rivermouth. The peak discharge of the flood exceeded the design discharge at that time and a large
bed scouring occurred in the river mouth. During the flood, temporal data of water levels
were measured at many observation points as shown in Figure 2. The cross-sectional bed
forms were surveyed before and after the flood.
14000

Discharge（m3/s）

12000
10000

Design discharge at that time
(9,000m3/s)

8000
6000
4000
2000
0
8/4 0:00

8/5 0:00

8/6 0:00

8/7 0:00

8/8 0:00

8/9 0:00

Figure 1: Discharge hydrograph of August 1981 flood of the Ishikari River mouth

Sea
Bed elevation(T.P.m)Flow
Observation point of water level

Figure 2: Plan-form and contour of bed level of the Ishikari River mouth and observation points of1981 flood
2.2 Previous Researches

There have been a lot of researches on the bed variation of the flood. One of them, Shimizu etal. (1986) developed 1D quasi-steady flow and bed variation analysis considering theconcentration of suspended sediment whose vertical distribution is assumed to be exponential
curve. They indicated that suspended load is dominant for the bed variation in the Ishikari
River mouth. Inoue et al. (2004) extended the above analysis to 2D quasi-steady flow and bed
variation analysis considering the effect of secondary flow. Where, suspended load is
calculated by concentration of suspended sediment and horizontal velocity whose vertical2

distributions are assumed to be exponential and polynomial curves, respectively. But a
problem to be solved remains that the bed scouring in meandering section was not explained
well because the secondary flow due to meandering channel are not solved well by 2D flow
analysis and because of the assumption of vertical distribution of suspended sedimentconcentration.
3.

The Analysis of Bed Variation during 1981 flood in the Ishikari River Mouth

3.1 Computational Method

We develop a new method of the unsteady numerical analysis by solving 1981 flood flow andbed variation so as to coincide with observed temporal changes in water surface profiles. The
discharge hydrograph and bed forms after the flood are compared with observed data. The
validity of calculated temporal bed elevations during the flood are verified by temporalchanges in observed water surface profiles. The analysis consists of quasi-3D unsteady flow
analysis (Uchida and Fukuoka, 2009) as shown in Eq. (1) ~ (6) and 2D bed variation analysis
for graded sediment (Fukuoka et al., 1998), employing the general coordinate system. The
bed load is calculated by the formula of Fukuoka (2010).
In the quasi-3D unsteady flow analysis, the horizontal velocity u i at   z s  z  h is assumed
as shown in Eq. (1). Where i,j = x,y, ui = the difference between depth averaged velocity U i
and bottom velocity ubi . ui is calculated by Eq. (2) using depth integrated vorticity  i .
ui  U i 

ui
2

2
3

1  3 ui  U i  ubi    ij 3  j h (1)2

Wh  2   ij 3  j hxi  3(2)

Where spatial change of depth averaged vertical velocity W is assumed to be negligible inshallow flows. Depth averaged horizontal velocity U is calculated by shallow waterequations as shown in Eq. (3) ~ (5) employing the general coordinate system (  , ).
J

h

(3)


 z
z 
U 
U 
U  ~   
 U  h ~  U  h ~  J U   U  cos   U  h ~  U  h ~    gh ~s  cos   ~s 
t

 
 



 



  0 hh U  h U  h


0
t

















z
z 
U 
U 
U  ~   
 U  h ~  U  h ~  J U   U  cos   U  h ~  U  h ~    gh cos   ~s  ~s 
t

 


  

  0 

(4)

1

~
  h~   J~h  ~ cos    ~ ~   ~ cos    ~ ~    h   J  

 









(5)



1

~
  h~   J~h  ~  ~ cos   ~   ~  ~ cos   ~    h   J  

 












 i is calculated by horizontal vorticity equation as shown in Eq. (6).
Jh i h Fi  Di  h Fi  Di  J ER zi  Pi 
t



3

(6)

Where Fi , Fi = flux of  i with expansion, contraction and rotation, Di , Di = flux of  iwith diffusion, ER zi = rotation term, Pi = production term.
3.2 Computational Method of Suspended LoadAt the section in question, effects of suspended load are dominant for bed variation andvertical distributions of concentrations of suspended sediment are affected by the longitudinal
deformation of channel shape and the secondary flow due to meandering channel. In thispaper, we compare results of two different methods for calculating suspended sedimenttransport.
At first (as calculation case 1), the vertical distributions of concentration of suspended
sediment is assumed as the exponential curves indicated by Lane-Kalinske (1950) as shown in
Eq. (7). Vertical distributions of horizontal velocities are assumed as polynomial curves asshown in Eq. (8) in the quasi-3D unsteady flow analysis. Depth integrated fluxes ofconcentrations of suspended sediment are shown in Eq. (9).
 w z
ck ( z )  cbk exp   0 k 
 


(7)

u i  3(8)

u i ( z)  U i 
u i ck 62 z  z  22  h2hcbk  2
1 h i
u i i exp        h   U i  u i 


u
c
dz



U


k


0
h
 
2 


(9)

Where z = height from reference plane, ck (z ) = concentration of suspended sediment ofdk(arbitrary grain diameter) at z, cbk = reference concentration of suspended sediment of dk,
w0 k = falling velocity of grain,  = coefficient of diffusion, u i (z ) = horizontal velocity( i   , ) at z, U i = depth averaged velocity, u i  U i  ub , h = water depth,i  3u i / h 2，   h /    s / w0k . Transport of suspended sediment concentration is
calculated by horizontal 2D advection-diffusion equations as shown in Eq. (10).
C k
C   
C 
1    u  c k   u  c k     k    k   ~  s ~  ~  s ~ 
t
J 




(10)

  
C   C    C 
C   cos    ~   s  ~k   ~   s  ~k   ~s  ~k ~  ~k ~   q suk  cbk w0 k
         J   

Where C k = Depth averaged concentration of suspended sediment of dk. Pick up rate ofsuspended sediment from bed surface ( q suk ) is calculated by the formula of Itakura and Kishi(1980).
Next (as calculation case 2), in order to calculate the suspended load affected by thesecondary flow directly, the concentration of suspended sediment is calculated by 3D
advection-diffusion equations as shown in Eq. (11). Where, the velocities are calculated byquasi-3D unsteady flow analysis. This analysis can calculate vertical distribution ofsuspended sediment concentration due to secondary flow.

4





c k 1    u  c k   u  c k  c k u    t  w0 k
c   
c 
 

 

 ~   s  ~k   ~   s  ~k 

t
J



 
 
   (11)  
c   c    c 
c    c  cos    ~   s  ~k   ~   s  ~k   ~s  ~k ~  ~k ~     s  k   q suk      z z     J   

3.3 Computational Conditions

The initial bed forms before 1981 flood were surveyed in 1979 and shown in Figure 2. The
boundary conditions of upstream and downstream ends are given by observed water level
hydrograph at 15km point and tide levels of the Otaru tidal observatory near the river mouth,
respectively. The mean grain diameter is 0.3mm downstream of 8km and 0.8mm upstream of
8km. The Manning’s roughness coefficients are determined so as to minimize the difference
between temporal changes in the observed and the calculated water surface profiles. As aresult, n=0.013 for the main channel downstream of 8km, n=0.021 for the main channel
upstream of 8km and n=0.050 for flood plains are given.
3.4 Computational Results

Water level (T.P.m)

Figure 3 shows temporal changes in observed and calculated water surface profiles in flood
rising and falling periods. Calculated water surface profiles coincide with observed waterlevels affected by bed variations in each time and there are little differences between
calculated water levels of each calculation cases. Figure 4 shows the comparison between
observed and calculated longitudinal bed forms before and after the flood. In the Ishikari
River mouth, suspended load causes the bed variations. The width-averaged bed elevations
after the flood of each calculation cases are similar and nearly coincide with that of observed
result. On the other hand, clear differences between each calculation cases can be seen in thelowest bed elevations of the meandering section from 3km to 8km.
8
7
6
5
4
3
2
1
0

left

Obs.
Cal. Case 1
Cal. Case 2

8/5 12:00
8/5 0:00-3k -2k -1k 0k(T.P.m)(T.P.m)levellevel
Water
Water

(Flood peak)
8/6 3:00

right

1k

2k

8
7
8
6
left right
7
5
Obs.
6
4
Cal. Case 1
5
Cal. Case 2
3
4
2
3
1
2
0
1
0 -3k -2k -1k 0k 1k 2k-3k -2k -1k 0k 1k 2k

3k

4k

5k

6k

7k

8k

9k 10k 11k 12k 13k 14k 15k 16k(a) Flood rising period(Flood peak)
8/6 3:00
8/6 18:00
8/7 12:00

3k

4k

5k

6k

7k

8k

9k 10k 11k 12k 13k 14k 15k 16k

3k

4k

5k

6k

7k

8k

9k 10k 11k 12k 13k 14k 15k 16k(b) Flood falling period
Figure 3: Temporal changes in observed and calculated water surface profiles

5

0
Averaged bed elevation

Bed elevation (T.P.m)

-2
-4
-6
-8

Lowest bed elevation

-10
-12

Averaged bed elevation-14

Lowest bed elevation

before flood Obs.
after flood Obs.
after flood Cal. Case1
after flood Cal. Case2

-16
-18

before flood Obs.
after flood Obs.
after flood Cal. Case1
after flood Cal. Case2-20-3k -2k -1k 0k

1k

2k

3k

4k

5k

6k

7k

8k

9k 10k 11k 12k 13k 14k 15k 16k

Figure 4: Average and lowest bed elevations before and after 1981 flood

Figure 5 shows cross-sectional bed forms of observation and calculations in meandering
section from 3.5km to 6.5km before and after the flood. The bed scourings are observed near
the left banks at 3.5km and 4.0km. On the other hand, the observed cross-sectional bed slopes
are almost flat at 5.5km and 6.0km. The cross-sectional bed slopes of calculation case 1(using Eq. (10)) are steeper than those of the observation at 3.5km and 4.0km and thescourings occur near right banks at 5.5km and 6.0km. On the other hand, the cross-sectional
bed slopes of calculation case 2 (using Eq. (11)) nearly coincide with those of the observation.
Figure 6 shows vertical distributions of calculated lateral velocities and concentrations ofsuspended sediment at 5.5km cross-section in flood rising period. The secondary flows due tomeandering channel can be seen at 5.5km cross-section as shown in Figure 6(a). The flows
near the bed transport suspended sediment near the bed. Therefore, there are small
concentrations at left bank and large at right bank in calculation case 2 (the concentrations ofsuspended sediment are calculated by 3D advection-diffusion equations) as shown in Figure
6(b).
10

3.5km

5

Bed elevation (m)Bed elevation (m)

10

0
-5
before flood Obs.
after flood Obs.
after flood Cal. Case1
after flood Cal. Case2

-10
-15
-20
-100

-15
-20

0

100

200

300

400

500 (m)

50

150

250

350

450

550

650 (m)600

700

800 (m)

10

4.0km

5

Bed elevation (m)Bed elevation (m)

0
-5

-10

10

0
-5

-10

6.0km

5
0
-5

-10

-15
-20
100

5.5km

5

-15

200

300

400

500

600

700

-20
200

(m)

300

400

500

Figure 5: Cross-sectional bed forms from 3.5km to 6.0km cross section before and after 1981 flood
10

10

5.0km

0
-5

0
-5

-10

6 -10

-15

-15

-20
100

200

300

400

6.5km

5
標高(m)

標高(m)

5

500

600

700

(m)

-20
300

400

500

600

700

800

900

(m)

4

Elevation (T.P.m)

2
0
-2
-4
-6
-8

water level
bed elevation
lateral velocity

velocity
0.2 m/s

-10
-12
100

150

200

250

300

350

400

(m)

(a) Vertical distributions of calculated lateral velocities
4
Elevation (T.P.m)

2
0
-2
-4
-6
-8

concentration (Cal. case1)

Concentration
0.001

-10
-12
100

150

concentration (Cal. case2)
200

250

300

350

400

(m)

(b) Vertical distributions of calculated concentrations of suspended sediment
Figure 6: Vertical distributions of calculated lateral velocities and concentrations of suspended
sediment at 5.5km cross-section in flood rising period

The discharge hydrograph is calculated by using temporal change of water surface profiles asshown in Figure 7. The accuracy of the observed discharge is not so high because large
temporal changes in bed elevations are not considered. Calculated discharge hydrographseems to be reliable more than observed one.
As above, the analysis evaluates properly the flood flow and the bed variation during 1981
flood of the Ishikari River mouth.
14,000Obs. (15km)

12,000Cal. Case1 (15km)
Cal. Case2 (15km)

Discharge (m3/s)

10,000
8,000
6,000
4,000
2,000
0
8/4 0:00

8/5 0:00

8/6 0:00

8/7 0:00

8/8 0:00

8/9 0:00

Figure 7: Discharge hydrograph at 15km upstream from the river mouth7

4.

Conclusions

In this paper, we develop the computation method of the quasi-3D unsteady flow and 2D bed
variation and calculate the August 1981 flood of the Ishikari River using observed temporalchanges in water surface profiles.
The bed variation mechanism during the 1981 flood can be evaluated by the present analysis
which explains well the discharge hydrograph during the flood period and bed elevations
surveyed after the flood.
Longitudinal changes in the cross-sectional bed forms at meandering section are calculated
appropriately by the present analysis. The main reason is because the suspended loadtransported by the secondary flow in the meandering section is solved well by the analysis of
the quasi-3D unsteady flow and the 3D advection-diffusion equations of suspended sedimentconcentration.
References
Fukuoka, S., Watanabe, A. and Okada, S. (1998) Bed topography analysis in a compound
meandering channel by using 3D numerical model with approximation of hydrostatics
pressure. Annual Journal of Hydraulic Engineering, Vol.42, 1015-1020, JSCE.
Fukuoka, S. (2010) River engineering adaptations against the global warming – Towards
generalization of close-to-nature rivers. Proceedings of the Japan Society of Civil
Engineers, F, Vol.66, No.4, 471-489, JSCE.
Fukuoka, S. (2011) What is the fundamentals of river design – Utilization of visible
techniques of sediment laden-flood flows. Advances in River Engineering, Vol.17, 8388, JSCE.
Inoue, T., Hamaki, M., Arai, N., Nakata, M., Takahashi, T., Hayashida, K., and Watanabe, Y.
(2004) Quasi-three-dimensional calculation of riverbed deformation during 1981-year
flood in the Ishikari River mouth. Advances in River Engineering, Vol.10, 101-106,
JSCE.
Itakura, T. and Kishi, T. (1980) Open channel flow with suspended sediments. Journal of
Hydraulics Division, Proceedings of ASCE, Vol.106, HY.8, 1325-1343.
Shimizu, Y., Itakura, T., Kishi, T. and Kuroki, M. (1986) Bed variations during the 1981
August flood in the lower Ishikari River. Annual Journal of Hydraulic Engineering,
Vol.30, 487-492, JSCE.
Takagi, J., Makino, N., Takemoto, N. and Morita, Y. (1982) Field survey of flood flow and
bed evolution in the lower Ishikari River. Annual Journal of Hydraulic Engineering,
Vol.26, 57-62, JSCE.
Uchida, T. and Fukuoka, S. (2009) A depth integrated model for 3D turbulence flows using
shallow water equations and horizontal vorticity equations. 33rd IAHR Congress,
Water Engineering for a Sustainable Environment, 1428-1435.

8