流域規模での洪水流出および土砂流出特性について
Runoff and Sediment Transport Characteristics on the Catchment Scale清水康行1・舛屋繁和2
Yasuyuki SHIMIZU and Shigekazu MASUYA
12北海道大学大学院工学研究科教授
北海道大学大学院工学研究科修士課程

要旨

河川管理を行う上で，流域規模での降雨－流出量を把握することはもちろんのこと，近年では，ダム等の河川構造物の維持管理，さらには水域における生態系保全等の観点から，流域規模での土砂輸送量の把握も求められている．このことから，降雨－流出量と土砂輸送量を流域全体で一括に取り扱うことが可能なモデルが，実際の河川管理を行う上で非常に有用であると考えられる．
Wongsaらは流域全体での流水，土砂輸送の追跡に加え，斜面域での降雨による土砂生産量の算定が可能な物理モデルとして，流域を土砂生産域としての斜面と土砂輸送域としての河道に分けた，流域スケール数値モデルを提案している．本
研究では，Wongsaらのモデルを北海道日高地方厚別川流域に適用し，同一流域における出水規模の差異によるモデル出力の違いを明らかにすると共に，日高地方
厚別川における2003年8月洪水を対象とし，土砂輸送量の把握，斜面域および河道域における，粒径区分別の土砂輸送特性について検討を行った．
《キーワード：洪水流出；土砂輸送，土砂侵食，流域規模》

１． はじめに我国の河川では，強い降雨時には流域斜面からの土砂流入や河道内での土砂輸送により，著しい河床変動が発生している事例が多い．そのため，河川管理を行う上で，流域規模での降雨－流出量を把握することはもちろんのこと，近年では，ダム等の河川構造物の維持管理，さらには水域における生態系保全等の観点から，流域規模での土砂輸送量の把握も求められている．このことから，降雨－流出量と土砂輸送量を流域全体で一括に取り扱うことが可能なモデルが，実際の河川管理を行う上で非常に有用であると考えられる．
流域規模での土砂輸送モデルは，これまでに橋本ら1)や，市川2)ら等によって研究が行われている．また，斜面の土砂流出モデルについては，裸地からの土砂流出に関する研究（村本ら3），村上ら4)）や，土石流堆積物の浸食による土砂生産に関する研究（金屋敷ら5））など様々な研究がなされている．また，Wongsaらは流域全体での流水，土砂輸送の追跡に加え，斜面域での降雨による土砂生産量の算定が可能な物理モデルとして，流域を土砂生産域としての斜面と土砂輸送域としての河道に分けた，流域スケール数値モデル（Wongsa and Shimizu6）等）を提案している．このモデルは，斜面上において物理モデルに基づいた土砂流出量を算出し，河道内において土砂輸送量を算定するモデルである．
Wongsaらは，このモデルを用いて，石狩川流域や沙流川流域における洪水ハイドログラフの再現や土砂輸送量の把握を試みている．
本研究では，Wongsaらのモデルを北海道日高地方厚別川流域に適用し，同一流域における出水規模の差異によるモデル出力の違いを明らかにすると共に，日高地方厚別川における2003年8月洪水を対象とし，土砂輸送量の把握，斜面域および河道域における，粒径区分別の土砂輸送特性について検討を行った．

２． 対象流域
本論文では，2003年8月に記録的な豪雨による被害を被った，日高地方厚別川を対象流域とした．厚別川は，源頭部を沙流川の支川である額平川に接し，支川である里平川・比宇川・元神部左の沢川・受乞川と合流し太平洋へ注ぐ，流域面積266.42km2 ，河道延長37.4kmの二級河川である．最上流域は1000mを越える急峻な山地に囲まれ，流域の大部分が原生林である．流域の土地利用形態は主に畑作であるが，下流域では稲作も見られる．対象流域の概要を図-1に示す．実際に数値計算を行う際には，図-2のように流域全体をいくつかの小流域に分割してモデル化し，250mメッシュ数値地図から，斜面域の流域面積及び斜面勾配，河道域の水路長及び水路勾配を算出した．また，実際の計算においては，河口までを計算対象とし，下流端の境界条件として海面水位を与えた．
①⑤沙流川
札幌

厚別川

②③新和 ⑦⑨④⑥比宇川⑧⑩⑪⑬新和比宇川

豊田⑫豊田

Pacific Coast

図-1

対象流域

図-2モデル化された対象流域

３． 計算モデルの概要および計算条件(1) 計算手順
本研究で用いる計算モデルは，以下に示す手順によって降雨－流出量の算定および土砂輸送量の算定を行う．
a) 斜面域における降雨－流出量の計算．
b) 斜面域における降雨および流水による土砂生産量の計算．
c) 斜面からの流入量を入力値とした，河道内における洪水追跡計算．
d) 斜面からの流入量を入力値とした，河道内における土砂輸送量の計算．
(2) 斜面域の支配方程式
斜面からの流出過程は，主に表面流と中間流に区分される．本研究では，計算方法が簡潔であること，洪水減衰期の再現性が良いことから，斜面上での流れの連続式および運動方程式として，次に示す表面流・中間流統合型Kinematic wave式を用いる．

∂hsl ∂q sl+= re∂t∂XkS slhsl

(1)


,0 < hsl < γDγq sl = kSα (hsl − γD )m + sl hslγ, hsl ≧ γD(2)ここで，hsl：斜面上の水深，qsl：斜面上の単位幅流量，re：当該斜面における有効降雨量，k：透水係数，Ssl：斜面勾配，γ：空隙率，D：中間流出層厚さ，α：Ssl1/2/nsl，nsl：斜面表面におけるマニングの粗度係数，m：5/3，X：斜面流下方向座標，t：時間である．なお，添え字slは斜面上であることを表している．
斜面の浸食は，雨滴が地表面に衝突することによる斜面の浸食，表面流による斜面の浸食，浸食された土砂の沈降の３つの過程から形成されているものと考えられる．これら３つの過程を考慮した土砂濃度連続式を(3)式に，斜面上の各粒径の構成比率連続式を(4)式に示す．

∂ (c sli hsl ) ∂ (c sli q sl )+= Dri + D fi − Ddi∂t∂X∂p1  ∂ (q Bsli )δ sli ++ Dri + D fi − Ddi  = 0∂t
1 − λ  ∂X


(3)
(4)

ここで，δ：交換層厚，λ：空隙率，csl：斜面流の土砂濃度，psl：斜面域の粒径iの構成比率，qBsl：斜面上の掃流砂量，Dr：降雨による生産土砂離脱率，Df：表面流による生産土砂離脱率，Dd：生産土砂沈降率である．なお，添え字iは混合粒径計算時におけるi番目の粒径に対する量であることを示している．なお，本研究では便宜的に，降雨期間中の山腹斜面の標高変化は無視している．
山地流域では，斜面上に樹木が生育している場合が非常に多い．そのため，裸地とは異なり樹木の遮断効果によって雨滴が直接地表面に到達する割合は小さく，降雨による生産土砂離脱率も低減されるものと考えられる．そこで本研究では，降雨による生産土砂離脱率として，次式（Wicksら9））を用いる．

Dri =k r Fw

ρs(1 − C )[(1 − C )M + M ]gcrd(5)

ここで，kr：降雨による生産土砂離脱係数，Fw：表面流水深補正係数，ρs：土粒子の密度，Cg：樹木の地表被覆率，Cc：樹冠の地表被覆率，Mr：降雨による単位面積あたりの運動量，Md：樹冠からの落下滴による単位面積あたりの運動量である．(5)式には，粒径毎に異なるパラメータは含まれていないことから，降雨による土砂離脱率は各粒径とも同じ計算値となる．しかし，(3)，(4)式において土砂沈降率Ddが考慮されているため，実際に流出する土砂量は各粒径毎に異なる．表面流による生産土砂離脱率
および生産土砂沈降率には，それぞれ(6)式，(7)式を用いる．

 k f  τ *csli  u*csli 1 − ,τ *sli > τ *csli 1 −(6)D fi =  ρ s  τ *sli u*sli ,τ *sli < τ *csli
0(7)
Ddi = k d w fi c sli
ここで，kf:表面流による生産土砂離脱係数，τ*sl：斜面上の表面流による無次元掃流力，τ*csl：斜面上の表面流による無次元限界掃流力，u*sl：斜面上表面流の摩擦速度，u*csl：斜面上表面流の限界摩擦速度， kd:生産土砂沈降係数， Wf:生産土砂沈降速度である．生産土砂沈降速度の算出には，次式に示す
Rubeyの沈降速度を用いる．
2w fi =

 6ν 26ν
sgd i +   −di3 di (8)ここで，s：砂の水中比重，d：粒径，ν：動粘性係数である．
(3) 河道域の支配方程式
河道内の流れに対する連続式および運動方程式は，そろぞれ(8)式，(9)式を用いる．

∂A ∂Q+= qL∂t ∂x∂Q ∂  Q 2  ∂H q Q
 + gA+Sf  = L 2+ ∂t ∂x  A  ∂x
 gA(9)(10)

ここで，A：河道の流水断面積，Q：河道内流量，qL：単位流路長当たりにおける斜面からの流入流量，
H：水位(H=η+h)，η：河床高，h：水深，Sf：摩擦勾配，g：重力加速度，x：河道流下方向座標である．
また，摩擦勾配はマニング則に従い，次式によって算出する．

 qn S f =  5/3 h 2(11)ここで，q：河道内単位幅流量，n：河道内におけるマニングの粗度係数である．
混合粒径計算時の河道内の土砂輸送を表す連続式は，粒径別浮遊砂濃度連続式と全流砂量連続式，各粒径の構成比率連続式が必要となる．本研究では，粒径別浮遊砂濃度の水深方向分布式として指数分布式を採用し，水深方向平均浮遊砂濃度の算出に(12)式を，粒径別浮遊砂濃度連続式として(13)式を用いる．

< ci >=cbi

β[1 − exp(− β )]∂(< ci > h ) 1 ∂(Q < ci > )q c+= qsui − w fi cbi + L bLi
∂t∂t
B
B

(12)
(13)

ここで，＜ci＞：水深方向平均浮遊砂濃度，cbi：基準点浮遊砂濃度，β：wfih/ε，ε：κu∗ιh/6，κ：カルマン定数，B：河道幅，qsui：河床からの土砂浮上量，cbL：単位流路長当たりにおける斜面からの流入浮遊砂基準濃度である．
混合粒径における各粒径の構成比率連続式は，交換層における粒径別の流砂の交換，河床変動および河床変動に伴って新たに交換層に加わる粒径分布を考慮し，次式で表される．

δ∂p iq c ∂η1  1 ∂ (q Bi B )+ pi*++ q sui − w fi cbi + L bLi  = 0∂t∂x 1 − λ  B ∂xB (14)ここで，δ：交換層厚，pi：河床材における粒径iの構成比率，λ：空隙率(0.4)，qB：掃流砂量である．
なお，河床変動の形態(堆積，掘削)によって，左辺第２項のpi*は異なる値を採用する必要がある．河床に土砂が堆積する場合は流砂中の粒度構成比率piを，河床に一旦堆積した土砂が掘削される場合は，土砂が堆積した時点における粒度構成比率pimを，元河床が掘削される場合には元河床の粒度構成比率pi0をそれぞれ採用する（清水7））．
全粒砂の連続式は，全粒径の流砂の交換を考え，次式によって表される．

∂η
1  1 ∂ ∑ i (q Bi B )∑ i (q L cbLi ) = 0++ ∑ i (q sui − w fi cbi ) +∂t 1 − λ  B
∂xB
(15)

ここで，Σiは全粒径の値を合計することを示す．
掃流砂量の算定には，芦田・道上の式8)を用いる．

q Bisgd i= pi 17τ *' i

3/ 2 τ *ci  u*ci 1 −1 −
u*i  τ *i (16)ここで，τ*'：河道内の無次元有効掃流力，τ*：河道内の無次元掃流力，τ*c：河道内の無次元限界掃流力，u*：河道内の摩擦速度，u*c：河道内の限界摩擦速度である．なお，混合粒径における各粒径iの限界摩擦速度は，次式に示す浅田の式9)によって算定する．
2
2u*ci log 23
 di=2u*cm   d i
  dm log 21 + 2    d m (17)ここで，添え字mは平均粒径に対する値であることを示す．
河床からの土砂浮上量算定には，次の板倉・岸10)の式を用いる．

 sgd i
qsui = pi K  α *Ωi − w fi u*′∞1ξ
exp − ξ 2 dξ∫′
τ *′i aτ′πΩi =+ *i − 1
B*i ∞ 1Bη
2
∫ exp − ξ dξ *i 0(a′π(18))()，

ξ *i =(19)ただし，

B*i = ξ i B*0

τ *′ciτ *′ci 0

ここで， τ*ci0 ：単一粒径で考えた場合の粒径 i の無次元限界掃流力， α* ： B*i/τ*'-1/η0 ， η0 ：0.5， K ：0.008，α*：0.14，Β*0：0.143である．
(4) 数値計算手法
河道部の流れの計算には，常流・射流が混在する流れでも再現性が良いCIP(Cubic InterpolatedPseudo-Particle)法を用い，土砂輸送連続の式には，それぞれ風上差分法を用いる．また，斜面部の流れ，斜面浸食の計算には風上差分法を用いる．
(5) 計算条件
本モデルでは，パラメータ以外の入力データとして時間降雨データが必要となる．本研究では各斜面にそれぞれ最も近い観測所の降雨データを入力時間降雨データとする．なお，入力されたデータは直接流出率による補正を行う．
各河道上流端では，境界条件として上流からの流入流量を与え，下流端では水位を与える．河道内の掃流砂・浮遊砂は，下流端から自由流出するものとしている．斜面域の境界条件については，各斜面上流端からの流入流量および流入流砂量は無いものとし，斜面下流端ではそれぞれ自由流出としている．また，計算初期条件は，河道内には基底流量が流れているものとし，斜面域の水深は0とした．
また，予備計算により，河道域でのパラメータはそれぞれ ∆x=500m， n=0.050～0.075，斜面域のパラメータはそれぞれ ∆X=100m， k=0.025， γ=0.20， D=0.4m， kr=20.0， kf=0.0005，Cg=0.8，Cc=0.8とした(Wicksら11），Lukeyら12））．なお，土砂粒径は最小0.01mmから最大400mmまでの粒径を15分割して土砂輸送の計算に用いた．

４． 計算結果(1) 洪水流量の再現結果
本研究では，前述した計算条件を用いて厚別川流域で発生した2003年8月洪水の解析を行った．その結果を図-3に示す．また，同じパラメータを用いて2001年9月洪水（図-4），2002年7月洪水（図-5）の解析を行った．なお，以降の図中添え字simは計算値を，添え字meaは実測値をそれぞれ示す．再現結果によると，2003年8月洪水および2001年9月洪水の再現計算では，洪水立ち上がり部および減衰部共に，洪水波形をある程度再現することができた．しかし，2002年7月洪水の再現計算においては，洪水波形の再現性が全く得られなかった．
2001年9月洪水および2003年8月洪水は，それぞれが記録的な豪雨であり，豊田観測所地点におけるピーク比流量が2003年8月洪水で9.4，2001年9月洪水では4.0であった．同じパラメータでの解析が可能であったことを鑑みると，両洪水の流出特性が酷似していたことが伺える．一方，2002年7月洪水は豊田観測所におけるピーク比流量が0.8程度であり，他の２洪水に比べて小規模な洪水であった．そのため，同じパラメータを用いての解析が行えなかったものと考えられる．
表-1 解析対象洪水の総雨量および流出係数
ピーク比流量（豊田観測総雨量(mm)所）
2001/9/10
246.8
4.0
2002/7/10
94.8
0.8
2003/8/8
327.7
9.4(2) 土砂輸送の再現結果
2003年8月洪水における，豊田観測所地点での土砂輸送量の再現結果を図-6に，豊田観測所地点における粒径別浮遊砂量を図-1に，流域全体における斜面域から河道域への流水および粒径別土砂流入量を図-8にそれぞれ示す．
図-6によると，掃流砂量，浮遊砂量の波形は，洪水流量の波形とほぼ一致していることがわかる．
また，図-8を見ると，斜面からの流水および土砂の流入量の波形が概ね一致していること，斜面域から河道域へ流入している土砂の大部分が，粒径0.075mm以下の細粒分であることがわかる．さらに，
図-7によると河道内の浮遊砂量は，そのほとんどが粒径0.075mm以下の細粒分である．以上のことから，河道域における浮遊砂量は，斜面域から河道域へ流入する細粒土砂から大きな影響を受けているものと考えられる．
出水前と出水後の河床材料粒径分布を比較すると（図-9），流域上流部の河道では河床材料が粗粒化し，下流部では細粒化していることがわかる．上流域では，斜面域から河道域への土砂流入量が比較的小さく，河床へ沈降する細粒土砂量よりも，河床から流出する細粒土砂量のほうが大きかった
001200105015006010007080500901000111213141time(hr)

51

6171

図-3 2003年8月洪水における洪水流量の計算値と実測値の比較（豊田観測所）

discharge(m3/sec)

30
40

3disharge(m /sec)

2000

10rain(m/sec)
Qmea(m3/sec)
Qsim(m3/sec)

1000

20rain(m/sec)rain(m/s)
Qmea(m3/s)
Qsim(m3/s)

2500

800

2030
4050

600

60
400

rainfall(m/sec)

3000

70
80

200

90
100

011121314151

61

71time(hr)

図-4 2001年9月洪水における洪水流量の計算値と実測値の比較（豊田観測所）

ものと考えられる．一方，下流域においては上流域から流出した細粒土砂に加え，斜面域から流入した細粒土砂が堆積し，河床材料の細粒化が発生したものと考えられる．また，土砂輸送計算による流域全体での土砂収支は，斜面域から河道域への土砂流入量が約2,180,000m3,河口から流出する浮遊砂量が約180,000m3，掃流砂量が約380,000m3となった．このことからも，下流域において細粒土砂の堆積が発生していることが推測される．

9846150810100

121450162011121314151

61

2500

76SS(m3/sec)
qb(mm3\sec)
Qsim(m3/sec)

541500
1000

3
2500
0
171

1121314151

61

71

time(hr)time(hr)

図-5 2002年7月洪水における洪水流量の計算値と実測値の比較（豊田観測所）

図-6 2003年8月洪水における土砂量計算値と流量計算値の比較（豊田観測所）

6.0

5.0

5.0

0.075mm以下
0.075mm～2.0mm
2.0mm以上
合計浮遊砂量4.0

slope erosion(m3/sec)

SS(m3/sec)

2000

10

180

3000

3.0
2.0
1.0

3500

4.0
3.00.075mm以下
0.075mm～2.0mm
2.0mm以上

3000

合計土砂流入量
runoff(m3/sec)

2000

2500

1500

2.0

1000
1.05000.00.011121314151

61

discharge(m3/sec)

102rainfall(m/sec)

discharge(m3/sec)

200

0

71

0
1

1121

314151

61

71

time(hr)time(hr)

図-7 2003年8月洪水における粒径別浮遊砂量（豊田観測所）

slope runoff(m3/sec)rain(mm/hr)
Qmea(m3/sec)
Qsim(m3/sec)

sediment
discharge(m3/sec)

250

図-8 2003年8月洪水における斜面域から河道への流水および粒径別土砂流入量（全河道区分）

図-9 出水前後における河床材料粒径分布の変化

５． 結論
本論文では，日高地方厚別川を対象流域とし，出水規模の差異によるモデル出力の違い，土砂輸送量の把握，斜面域および河道域における，粒径区分別の土砂輸送特性について検討を行った．以下に，本論文で得られた結論を列挙する．
a) 流出特性が類似していると考えられる洪水では，同じモデルパラメータを用いた流出解析がある程度可能であるが，それ以外の洪水を解析する場合には，新たにモデルパラメータの検討が必要とな
る．
b) 浮遊砂として河道に存在している土砂および斜面域から河道域へ流入する土砂は，その大部分が粒径0.075mm以下の細粒分である．
c) 斜面域から河道域へ流入する細粒土砂は，河道内の浮遊砂量に対して大きな影響を与えている．
今後，より多くの流域に本モデルを適用することにより，モデルパラメータの決定方法を確立すると共に，さらに詳細な土砂流出特性の検討が必要であると考えられる．

参考文献
1)

橋本晴行・朴埼璨・高岡広樹・荒渡光貴：山地流域における豪雨時の水と土砂の流出解析，水工学
論文集，第 47 巻，pp.745-750，2003.

2)市川温・藤原一樹・中川勝広・椎葉充晴・池淵周一：沖縄地方における赤土流出モデルの開発，水工学論文集，第 47 巻，pp.751-756，2003.

3)

村本嘉雄・道上正規・下島栄一：大戸川における微細砂の流送過程について，京都大学防災研究所
年報，No.16B，pp.433-447，1973.

4)村上正吾・林誠二・亀山哲・渡辺正孝：河川流域内農林地域における土砂動態モデルに関する基礎的検討，水工学論文集，第 45 巻，pp.799-804，2001.

5)

金屋敷忠義・芦田和男・江頭進治：山地流域における濁度物質の生産・流出モデルに関する研究，
第 24 回水理講演会論文集，pp.143-151，1980.

6)

Wongsa, S. and Shimizu, Y.：Modelling Pre-Channelization and their Impact on Flood and Sediment Yield inIshikari River Basin, Ann. J. of Hydraulic Eng, Vol.47, pp.223-228, 2003.

7)清水康行：沖積河川の縦断形と河床材料分布形の形成について，土木学会論文集，No.521/II-32，pp.69-78，1995.

8)

芦田和男・道上正規：混合砂れきの流砂量と河床変動に関する研究，京都大学防災研究所年報，
No.14B，1971.

9)

浅田弘：山地河川の流砂量と貯砂池の堆砂過程に関する研究，電力中央研究所報告総合報告，2，1976.

10) Itakura, T. and Kishi, T.：Open Channel Flow with Suspended Sediments, Proc. of ASCE, No.106(HY8),pp.1325-1343, 1980.
11) Wicks, J. M., and Bathurst, J.C.：SHESED: A Physical Based, Distributed Erosion and Sediment Yield
Component for the SHE Hydrological Modeling System. J. Hydrol., 175, pp.213-238, 1996.
12) Lukey, B.T., Sheffield, J., Bathurst, J.C., Hiley, R.A. and Mathys, N.：Test of the SHETRAN Technology forModeling the Impact of Reforestation on Badlands Runoff and Sediment Yield at Draix, France, J. of Hydrol.,235, pp.44-62, 2000.