こうえいフォーラム第 13 号 / 2005. 1

河口土砂流検討に平面2次元移動床解析モデル-NKhydro2Dの適用
APPLICATION OF DEPTH-INTEGRATED MODEL-NKhydro2D TO THE EVALUATION OF BED
DEFORMATION IN RIVER ESTUARY

金海生*・千田健一**・佐々木成人**・山田裕康***

Haisheng JIN,Kenichi CHIDA, Shigeto SASAKI and Hiroyasu YAMADA

A method using the NKhydro2D model for predicting flow and sediment transport in river estuaries isreported in the paper. The model solves the depth-integrated 2D Reynolds' equations which are based onorthogonal curvilinear coordinate system. Both bed load and suspended load transport are taken into accountin bed deformation. The method was applied to simulate the flow and sediment transport in a river estuary.
The computational results of water surface level, flow velocity, bed deformation, and sediment transportvolume are compared with the measured data during a flood. It shows that the method predicted very well theflow and sediment transport in the river estuary. It is anticipated that the method can be applied to issues
such as sediment control and training works in river estuaries.
Key Words： Bed deformation, bed load, suspended sediment, sediment transport volume, velocity,
river estuary, NKhydro2D model, boundary-fitted coordinates

1. はじめに
国土資産の大半が河川下流の沖積平野に集中している日本では、河口部における洪水疎通能力の確保が治水上極めて重要である。一般的に、浚渫などを行っても、時間の経過とともに土砂堆積が生じ、洪水流下能力を確保できない場合が多い。安全性の高い河口計画を行うためには河口部地形変化の予測技術の確立が強く求められている。
河川流が卓越する河口部の河床変動に関しては、近年、理論的・実験的な研究により、浮遊砂が支配的な場合、流軸では侵食、両側に自然堤防のような堆積が生じ、鳥趾状デルタに近い堆積形状が発達することが示されている。浮遊砂の濃度と、流速の流下方向への逓減の仕方によって河口から遠い領域での流軸上の堆積形状が決定される。浮遊砂の移流拡散の影響が小さく浮遊砂濃度の逓減が流速より速いとき、河口から遠い領域でも侵食傾向となるが、流速の逓減が浮遊砂濃度より速いとき、河口から遠い領域では流軸上に浮遊砂の堆積が生じることになる。一方、掃流砂が優勢のとき、流軸上に堆積そしてその両側に侵食が生じ、円弧状のデルタが形成される1）。

口工事、例えば、導流堤の建設などによりこのバランスが崩れ、堆積が優勢となった場合、河床上昇により洪水流下断面が不足し、計画高水以下の流量を安全に流下できず、出水時に河川水位が計画高水位を超える可能性がある。逆に、洪水時河床低下により導流堤付近の河床が異常に侵食された場合は、導流堤の構造物としての安定性が失われる可能性がある。
留萌川は、留萌港内に流入しており、河口から流出した土砂が港内に堆積するため、留萌港内の泊地で航路の水深確保を目的として維持浚渫が行われているが、導流堤の建設（移設）により、港内への留萌川河口流出土砂の量・堆積分布が変化すると予想される。この変化を把握する必要がある。
本研究では、境界適合直交曲線座標系に基づいた水深積分の平面2次元移動床水理解析モデル−NKhydro2Dモデル2）により留萌川河口部のフラッシュシミュレーションを実施
し、洪水期河口の流れ、掃流・浮遊土砂輸送および河床変動の評価を行った。ここに、実績洪水データを用いたモデルの検証結果を報告する。

したがって、河口付近の河床高・河道断面形は、河口砂
州の延伸、漂砂、河川上流から運ばれた土砂の堆積を含む

2. 水深平均の2次元移動床水理解析モデル

堆積過程と、常時の流れによる侵食、洪水時河口フラッシ

水理解析においてはできるだけ実河川に近い河道条件（平ュを含む侵食過程、のバランスの上に成り立っている。河面形状、抵抗など）を与えるのがベストである。そこで、境界適合型の直交曲線メッシュをベースにした移動床水理解析* 中央研究所 総合技術開発部** 札幌支店 技術部*** 福岡支店 技術部

モデル−NKhydro2Dを適用し、流れおよび流送土砂の量・侵食 / 堆積分布がどのように変化するかを解明することとした。

55

河口土砂流検討に平面2次元移動床解析モデル-NKhydro2Dの適用（1） 直交曲線座標系における流れの基礎式：渦（有効）粘性係数

水深平均のReynolds方程式を直交曲線座標系に変：摩擦速度

換すると（1）〜（3）式になる。流れの支配方程式には植生などの影響を考慮する。高木などによる抗力の影響につ：底面せん断力

いては運動量保存式（2）〜（3）の：方向の底面せん断力：方向の底面せん断力項を参照されたい。

その他の植生の影響は植生の種類、密生度などの特性に応じてマニング粗度係数に考慮する。
＜水の質量連続式＞
（1）
：水の密度

＜流れの運動量方程式＞：重力加速度：マニング粗度係数（河道計画によって河床材料から決まる低水路粗度係数、植生などから決まる高水敷の粗度係数（2）を基準に設定するか実測流速・水位データなどを用いてキャリブレーションする）：植生（高木などによる）の抗力係数：植生による単位体積あたりの遮蔽面積：直交曲線座標：デカルト座標：座標変換における係数5）（3）

ここに、：方向の流速：方向の流速：デカルト座標系における：デカルト座標系における（2）河床変動に関する支配方程式
河床変動は流砂輸送によるものである。河床高の時間的な変化は次の流砂の連続式を用いて表される。
方向の流速
方向の流速

：水深
：水位および河床高（解析格子点におけ（4）
ここに、は河床高、径、形状などに依存する）、は河床材料の空隙率（砂の粒はそれぞれ方向のる河床高は、対象区間の横断測量成果単位幅、単位時間あたりの体積流砂量である。本検討では、を内挿して与える）

砂の浮遊輸送も考慮し、流砂量を掃流砂輸送量： 水深平均Reynolds応力4）

輸送量と浮遊砂からなることとする。混合砂の場合、それぞれの粒径階の流砂量を求め、その合計を流砂量とする。
掃流砂輸送量およびその輸送方向は、局所的な河床勾配および流線湾曲による二次流の影響を考慮し、芦田・道上式に準じて求める 6）。この流砂量式は有効掃流力を考慮して掃流砂量を推算することとしている。
浮遊砂輸送量は浮遊砂濃度と流速の水深積分からなる。図−1は、流速および浮遊砂濃度の鉛直方向分布の概念図である。浮遊砂の場合、河床変動は沈降量と巻き上げ56

こうえいフォーラム第 13 号 / 2005. 1

量の差によって発生する。流れが速く巻き上げが優勢のと単位面積から単位時間に巻き上げる粒子群の体積き、河床侵食が生じ、逆に沈降が支配的なとき、堆積が生位：じる。巻き上げられた浮遊砂は移流と拡散によって横断方する粒子群の体積に等しいため、向に広がる。水深平均の浮遊砂濃度は式（5）の浮遊砂質量立つ。

での理論的な平衡濃度であり、経験保存式（移流拡散式）で解析することができる。

式によって計算することとなる。例えば、本研究では、板（単）は、単位面積当たり単位時間に河床へ沈降は基準点

倉・岸の式により計算する8）。なお、係数河床近傍の実際の浮遊砂濃度

ここで、ればが成りを0.0018とした。

は下記式で求める。

は水深平均の拡散係数であり、Prandtl理論によである。

：摩擦速度、：von Karman定数

単位面積あたり単位時間の砂の浮遊運動による河床変動計算においては以下の式が成り立つ。
砂の浮遊運動による＝（沈降による河床上昇）−（浮上に図−1 流速および浮遊砂濃度の鉛直方向分布の概念図よる河床低下）=
（5）

（7）

したがって、掃流砂とともに浮遊輸送が存在する場合の河床変動に関する支配方程式は式（8）に変更することができる。
（8）

ここに、度、は粒径は水深平均、は河床近傍の浮遊砂体積濃は河床面から砂粒子の巻き上げによる浮上量、の砂の沈降速度、は拡散係数、は乱流拡散Schmidt/Prandtl数である。混合砂の場合、それぞれの粒径階に対して、式（5）を用いて、浮遊砂濃度を解析する。

はvan Rijnの式7）を用いる。

上記は均一砂をベースに導いた式ではあるが、混合砂の場合も似たような方法をとることが一般的である。
混合砂礫を対象とする場合は、交換層における粒径階の含有率を当粒径階の砂礫の質量保存則により求める。

（3）初期条件および境界条件
解析開始時、各計算点における水位（水深）、流速および（6）

浮遊砂濃度（通常、初期0とすることが可能）を仮定し、下流端において計算開始時の水位、上流端において計算開始時の流量、浮遊砂濃度を与えて解析を実施し、各計算点の水位、流速および浮遊砂濃度が一定値に収束した状態を初

ここで、、は砂の密度、は水の動粘性係数である。

期条件とする。また、各計算点において、初期河床高および粒度分布を初期条件として与える。

単位幅、単位時間あたりの浮遊砂輸送量

は浮遊砂濃度流速の水深積分からなる。

境界条件としては、上流端において流量時系列および浮遊砂濃度、下流端において水位時系列を与える。境界での混合粒径ごとの掃流砂量は掃流力の計算値より求めた平衡

ここに、は浮遊砂濃度の基準点とされる河床面からの高さであり、の取り方により濃度分布が若干異なってくる。その結果、浮遊砂量の算定は、濃度基準点の点の濃度は水深の値とその算定問題に帰着する。一般的に、の5％河床面から砂の浮上量の値が採用されている。の計算としては、いくつかの経験式が提案されている。河床が平衡状態にあるとき、河床の掃流砂量を粒度分布にあわせて与える。
（4）解析スキーム
流速と水位、河床高、浮遊砂濃度変数をスタッガード格子上に配列し、有限体積法に基づいて解析を行う。空間微分項の離散では、移流項、レイノルズ応力項にPower Law
スキーム9）を、その他の項に中心差分を適用する。時間微57

河口土砂流検討に平面2次元移動床解析モデル-NKhydro2Dの適用

分項には、陰形式のスキームを用いる。

（1） 水位の検証

計算の手順としては以下のとおりである。

河口観測所断面における水位の計算結果（実線）と実測① 運動量式（2）、（3）と水の連続式（1）を結合することにより水位を求める。

値（三角）を図−5に示す。計算の流量（破線）も示している。

② 水位の計算結果をもとに運動量式（2）、（3）により流速の解析を実施する。

河口観測所断面において、洪水ピーク時の実測水位は6.0m、最大水深は約7.33mであった。本再現計算では実際③ 上記水理量をベースに浮遊砂の巻上げ量を求める。に発生した下流市街地の氾濫が考慮されていないため、計④ 浮遊砂の移流拡散式（5）により浮遊砂濃度を計算する。

算のピーク水位は若干高くなっている。しかし、計算のピ⑤ 水理量をベースに掃流砂量を算出する。

ーク水位は6.09m、最大水深は7.42mで、最大水深の相対誤⑥ 河床材料の連続式（4）または（8）により河床高変化差は1.2％（=水位差 / 水深×100％）である。

を求める。
⑦ 各粒径階

一方、中低水流量時の計算水位は実測値より若干低くなの質量保存則によりその含有率の変化を計算する。
以上の手順を繰り返して解析を進める。っている。原因は以下のように推測される。
a）用いた河道の河床地形データは昭和63年洪水後の測量データに基づいたため、洪水期河床侵食の影響がまだ残っており、その影響が中低水流量時のほうで大きか3. 留萌川河口への適用った。

前記平面2次元移動床解析モデル−NKhydro2Dを用いて、

b）実際の河道では中低水流量時に河床波が発達していた

留萌川の河道を含む河口域流れおよび河床変動のシミュレ可能性があるため、河床波による付加抵抗が生じ、今ーションを実施した。ここに、近年最大洪水の実績データ回の解析で見積もった抵抗（洪水中に一定とした）よに基いたモデル検証結果を報告する。

り大きくなり、水位が計算値より若干高くなった。

昭和63年8月、留萌川流域において大きな出水があった。

c）計算では河口観測所流量を使用している。河口観測所

KP2.8に位置する留萌河口観測所におけるピーク流量はが河口から約3kmに位置するため、それより下流区間880m3/s、ピーク水位は6.0mであった。この洪水のピーク後の流入が河川の水位に影響していた。その影響は、堤4時間ごろ、航空撮影により流速観測が行われた 10）。そこ

内地盤高より低い水位の中規模流量時に比較的に大きで、この洪水の観測データおよび既往調査の流砂量結果にかった。一方、洪水ピーク前後では、河川水位が高く、基づき、解析モデルの再現性を検証する。

区間の流入が少なく、また、一部の地域では氾濫が起

解析領域は、河道部KP5.0の大和田橋を上流端とし、下流は海域の防波堤範囲とした。河岸境界および防波堤など構こったため、計算水位が観測値より約0.1mほど高くなった。

造物に合った境界適合型直交曲線メッシュの一部を図−2d）減水時における河口観測所流量（急激に下がる）と水位に示す。メッシュサイズは横断方向に最小約2m、縦断方向（緩やかに下がる）の変化特徴から、以下のことが考えに最小約5mで、総メッシュ数は約10,000である。

られる。①洪水ピーク後、下流市街地に氾濫で溢れた水検証計算用の河床高および河床材料について、洪水前のの一部が河川へ還流したことにより、下流河川の水位が実測データが必要となるが、このようなデータは存在して上昇し、河口観測所断面の水位に影響したこと。②氾濫いなかった。したがって、昭和63年洪水に最も近い時期にを無視して、洪水が河川に集中して流れる条件で解析し調査したデータを用いることとした。すなわち、河床高は、たため、洪水ピーク前後の河床侵食が若干大きめに見積

河道部に昭和63年12月横断測量データ（洪水後）、海域部にもられ、洪水後の計算水位が低くなった。

昭和58年8月留萌港深浅測量図を用い、欠測部分を平成6年9以上の考察から、初期河床、河床波に伴う粗度評価、残月深浅図と昭和61年3月留萌港港湾計画図で補完した。河床

流域の氾濫を含む流出入などの条件をより正確に計算に反材料は、河道部に平成8年度河床材料調査データ 11）、海域映できたならば、計算水位は、中低水流量時も含めて実績部に平成6年度漂砂調査データ（図−3）を用い、11分級の水位により良く整合させることができたと考える。すなわ

混合砂として扱った。河道の下流および河口部の河床材料

ち、NKhydro2Dモデルは留萌川の実績洪水水位を十分に再に細粒分が多いことが判る。

現し得たと判断する。

洪水当時の河口観測所断面の流量ハイドロと留萌港の潮位観測データ（図−4）をそれぞれ上、下流の流れの境界条件とした。流砂量条件については、上流端平衡給砂とした。

58

こうえいフォーラム第 13 号 / 2005. 1

図−2 境界適合直交曲線メッシュの一部（2）流速の検証
図−6に航空写真撮影で測った洪水ピーク後4時間ごろの平均流速ベクトルおよび流速コンターを示す。また、洪水ピークから約4時間後の河口部航空写真を図−7に示す。
図−7の赤い矢印は洪水の主流部イメージであり、西と外北防波堤間の主流が西寄りになっていたことを示しているが、東寄りの分岐も無視できないようである。
航空写真調査により、河口から留萌港範囲の洪水流（流速、流向）について下記のような特徴が判明した。
a）流向は、河口から西防波堤、外北防波堤の間に向う流れが主流である。
b）導流堤を越流した流れは、古丹浜防波堤の開口部を経

59

河口土砂流検討に平面2次元移動床解析モデル-NKhydro2Dの適用て外港に至る。
c）南防波堤と北防波堤の間になる外港出口は複雑な流れとなっている。

響によるものと考えられる。
図−8、9に洪水ピーク4時間後の流速ベクトルおよび流速コンターの計算結果を示す。図−8、9に示す通り、導流d）流速は、河道部で約5m/sだが導流堤部で約3m/sに減少

堤部で河道流の主流は導流堤寄りに流れ、洪水ピーク時導している（横越流の影響）。河口部では、水面の急降下流堤に発生した横越流は、古丹浜防波堤の開口部を経て外とともに流れが加速流となり、流速が約4m/sに加速し、港に至り、河口から流れる主流と合流する。河道部で約流下とともに減速して沖合700〜800mで1m/sとなる。

5m/sあった流速は導流堤部で約3m/sに減少し（横越流の影e）導流堤部で河道流の主流は導流堤寄りに流れている。

響）、一方、河口部で水面の急降下とともに流れが加速流とf）港内の流向は、風向、波向により変化するが、大まかなり、出口直後に流速が約4m/sに達する。その後、主流はにみると右回りの環流を形成している。これは、南防減速しながら、河口から西防波堤、北外防波堤の間に向う。

波堤と西防波堤の開口部からの流入と港内の地形の影

港内の流向は、右回りの環流を呈している。これらの特徴
は図−6、7の洪水時航空写真撮影結果（上記a〜fにまとめている）とよく一致している。

（3）河床変動および浮遊砂濃度の検証
図−10に、洪水期間中の河床高変動量計算結果を示す。
河床変動は河口幅の数倍程度の範囲に集中的に起こっていることが判る。洪水期間中、河口部においては河床が約1m程度、局所的に2mほど（河口KP-0.6〜KP-0.4前後の左岸側導流堤沿いに）侵食された。河口フロントでは、洪水流心部に1mほどの侵食、その両側に0.5mほどの堆積が生じた。
上流から流送された土砂および河口から洗掘された土砂が河口沖合に最大約4m堆積した。留萌港内その他領域には堆
図−3 河道0〜5K区間および河口部の河床材料積が小さかった。この結果は定性的に泉ら 1） の研究成果

図−4 洪水流量ハイドロ、留萌港潮位（昭和63年8月）

図−5 河口観測所断面の水位の計算結果と実測値

60

図−6 昭和63年8月洪水流速ベクトルおよびコンター（航空写真撮影による、ピーク約4時間後）

こうえいフォーラム第 13 号 / 2005. 1（文献1の図−11、12）と一致している。
図−11に、昭和58年と平成6年の深浅測量データを基に作成した河床高変化のコンターを示す。その間、昭和63年8月に大出水があった。図−10、11からは、流軸での侵食、流軸横での堆積および河口沖合での堆積の解析結果と実測結果が定性的によく一致していることが判る。
図−12に洪水ピーク4時間後の浮遊砂SS濃度の計算結果を示す。図−7に貯木場と河口フロントからの赤い矢印に沿った明るい部分は濁り洪水流の主流部（すなわち土砂濃度の高い領域）であったことが推測される。図−12は上記範囲でSS濃度の高いことを示している。したがって、図−7と比較すれば、SS濃度分布の計算結果からもNKhydro2Dモデルは昭和63年8月洪水を良好に再現していることが判る。

図−8 洪水ピーク4時間後河口部流速ベクトルの計算結果

図−7 洪水ピーク4時間後の航空写真（赤い矢印は主流部）
図−9 洪水ピーク4時間後河口部流速コンターの計算結果
61

河口土砂流検討に平面2次元移動床解析モデル-NKhydro2Dの適用

図−10 洪水期間中の河床高変動量の計算結果（マイナスコンター＝河床侵食、プラスコンター＝河床堆積）図−12 洪水ピーク4時間後浮遊砂SS濃度の計算結果
砂輸送には以下の特徴がある。
a）全流砂量の内、浮遊砂量が掃流砂量より圧倒的に多い。
b）河川整備前の年平均流砂量（河口から南防波堤の北側末端に至る範囲の堆砂量ベース変換値）実績は約59千m3であった。
図−13、14および表−1によれば、空隙率を仮に0.4とした場合、洪水期間中に河口観測所断面を通過する土砂総量の計算結果は約145.8千m3であり、そのうち浮遊砂が約98％占めて圧倒的に多い。河口断面の通過土砂量は208千m 3であるが、そのうち、約半分以上の106.8千m3が粒径0.04mm以下の極端に細かい粘土質のもので、河口部に堆積の多いシルトより粗い砂分が約101.2千m3である。
粒径の細かいシルト質の土砂は一般に浮遊を続けながら港外を目指して拡散する。図−12のSS濃度分布の計算結果

図−11 昭和58年から平成6年までの河床高変化の測量結果（平成6年河床高−昭和58年。マイナス＝河床侵食、プラス＝堆積）（4）流砂量の検証
対象洪水期間中に河口観測所断面（KP2.8）と河口断面（KP-0.2）を通過する土砂量の計算結果をそれぞれ図−13、14に示す。表−1はその数値である。
一方、昭和50年代からの現地調査によれば、留萌川の流62がそのような拡散を再現している。実測堆砂量調査領域内の堆砂量の解析結果（堆砂量ベース）は109千m3である。
計算の港内堆砂量（一洪水に対する）109千m 3を直接比較する現地データはないが、実測の年平均堆砂量59千m3と比較した時、計算対象の洪水規模からして計算値が約2倍の量となったことは十分に納得できる。

こうえいフォーラム第 13 号 / 2005. 1d）浮遊砂濃度および河口部河床変動の現象が合理的に再現され、深浅測量結果および他の研究成果と一致した。
e）流送土砂の構成および流送土砂量の解析結果が留萌川
河口部土砂堆積量の調査結果と一致した。
河口工学は、従来半経験的に構築されてきたが、今回適用したNKhydro2Dモデルにより、河口地形の予測が可能となった。今後、当手法は以下の分野への応用が期待される。①河道、河口から沿岸にかけての土砂管理、②河口洪水流下断面の動的管理、③河口改修、導流堤レイアウト検討など。
図−13 洪水期間中河口観測所断面を通過する土砂量（空隙率無視。縦軸：左側＝掃流砂、浮遊砂の粗い分、右側＝細かい浮遊砂、シルト）

謝辞： 本研究の一部は、北海道開発局留萌開発建設部から受注した「留萌川導流堤施設調査検討業務」で実施されたものである。ここに、本成果とりまとめにあたり貴重な御助言と御指導をいただいた留萌開発建設部治水課の関係の方々に深謝する次第である。

参考文献
1）泉典洋、田中仁、坪井宏介、伊達政直：河川流が卓越する河口テラスの初期堆積形状に関する実験、土木学会論文集、No.740/Ⅱ-64、2003.
2）金 海生、佐々木成人、山田裕康：急流扇状地河川の河床変動および河岸侵食のシミュレーション、こうえいフォーラム、

図−14 洪水期間中河口断面を通過する土砂量（空隙率無視。縦軸：左側＝掃流砂、浮遊砂の粗い分、右側＝細かい浮遊砂、シルト）

No.12、2004.
3）Kalkwijk, J.P.Th., De Vriend, H.J.: Computation of the flow inshallow river bends, Journal of Hydraulic Research, Vol.18, No.4,

表−1

洪水期間中の通過土砂量（単位：1,000m3）

1980.
4）Rodi, W.: Turbulence models and their application in hydraulics,
IAHR Publication, DELFT, The Netherlands, 1980.
5）Thompson, Joe F., Warsi, Z.U.A., Mastin, C.W.:: Numerical Grid
Generation ― foundations and applications, Elsevier Science
Publishing Co., Inc., New York, 1985.

4. おわりに境界適合直交曲線座標系に基づいた水深積分の平面2次元移動床水理解析モデル−NKhydro2Dモデルを留萌川河道および河口部に適用した。近年最大洪水の実績データを対象にしてモデルの検証を行った結果、NKhydro2Dモデルは留萌川河口の流れおよび土砂移動を良好に再現することが判り、以下の結論を得た。
a）洪水期間中の河口観測所断面で水位の時間変化および洪水ピーク前後に導流堤部での急な水面勾配変化を高精度に再現した。
b）河道内導流堤部の減速流、河口部の加速流、導流堤の

6）芦田和男、道上正規：混合砂礫床の河床変動に関する研究、 京
大防災研究所年報、 No.14（B）、1971.
7）Van Rijn, L.C.: Mathematical modelling of suspended sediment in
non-uniform flow, Journal of Hydraulic Engineering, ASCE, Vol.110,
No.11, 1984.
8）土木学会：平成11年度版水理公式集、p.167、1999年11月9）Patankar, S. V.: Numerical Heat Transfer and Fluid Flow, Hemisphere
Publishing Corporation, 1980.
10）留萌開発建設部：昭和63年度留萌川・洪水測定業務報告書、
1989.111）留萌開発建設部、株式会社開発工営社：平成8年度留萌川改修工事の内留萌川河床材料調査業務報告書、1996.10

横越流などを合理的に再現した。
c）洪水流の主流部位置およびその流向、留萌港内全域に発生する右回りの環流などを良好に再現した。

63