土石流の Fortran コード

Hyper KANAKO の不具合が多く、代替えコードを使用することになりました。

先輩が調達した土石流の Fortran コードを読んでいたのですが、これがとても読み易い。素晴らしい書き方だと感心しながら読んでいました。

中身は高橋モデルの一部が使用されており、少し手を加えることで河川砂防技術基準に適合させることができそうでした。流れがシンプルかつ明快で、偏微分の離散化もわかりやすく、すぐに読み終わりました。

おそらく、KANAKO の計算部分も同様にシンプルなはずです。なのになぜ不具合が多いのでしょう？

おそらく後付けの GUI 側の不具合がメインなのでしょう。GUI はソフトウェアメーカーさんに任せるのが一番良いのですが、残念ながら FORUM8 さんの開発は止まっています。ユーザーが少なく採算に合わないので Ver.UP されないのだと想像されます。

エンジンと GUI を完全に切り離し、データフォーマットをマニュアルに載せておけば扱いやすくなるでしょう。

ま、Fortran コードも整備しながら、ゆっくり待ちましょう。

「高橋のモデル」平衡土砂濃度

平衡土砂濃度を３区分する水面勾配の閾値（0.03、0.138）と、平衡土砂濃度算出式の解釈に参考となる説明が以下にありました。※式2

高橋「土石流発生のメカニズム (2)」

https://www.jstage.jst.go.jp/article/sabo1973/45/3/45_3_38/_article/-char/ja/

式 (2) では, 勾配が急になればなるほど平衡砂礫濃度が大きくなるが, 密に詰めた充填濃度になれば, 実際上流動は不可能となる。実験によれば, 流動が可能な最高濃度は0.7程度で, 式 (2) がそれ以上の値となる時でも0.7あるいはもっと小さい濃度で流れる1)。一方, 式 (2) は,勾配が小さくなればいくらでも小さい砂礫濃度を与える。しかしながら, 実験によれば, 砂礫濃度は0.2～0.3以下では, 砂礫は全流動層には分散できず, 図-1 (2) に示しているような掃流状集合流動となる。図-2は通常の粗い砂礫を用いた実験によって得られた砂礫濃度と勾配の関係と, 式 (2), 掃流状集合流動に関して著者が理論的に求めた平衡濃度csooを与える式(3)2), および水山の経験式3), さらに, 芦田・高橋・水山の掃流砂量式4)による掃流砂の土砂濃度とを示している。砂礫濃度がいかに勾配に支配されて変化するかがよく分かる。図-2によれば,土石流と掃流状集合流動の境界は流路勾配10. 程度にあり, 十分長い一様な堆積層があれば, 100程度の渓床で土石流が発生できることを示している。なお, 式 (3) は以下のようである。

2013年 Hyper KANAKO 勉強会の第一部で紹介されたスライドによれば、掃流状集合流動と掃流砂の境界：0.03と土石流と掃流状集合流動の境界：0.138は、勾配1.7度、7.9度に相当します。図では3度、10度程度となっていますので、小さめの値が採用されています。その後の実験で決めた閾値だと推察されるのですが、出典がないため追いかけられません。

平衡土砂濃度（掃流状集合流動（式3）、掃流砂）の導出や出典はよくわかりませんでした。理論式だそうですが。

いずれにしても、ある理論式をある範囲で使い分けている、範囲の閾値は実験で決めているというようなことだと推察されます。

ひとまずここまで。

「高橋のモデル」侵食/堆積速度

まずは、侵食/堆積速度 i の確認です（特に侵食・堆積速度係数の一般値の根拠）。

高橋ほか「変勾配流路における土石流の形成」では、侵食速度は侵食深さ aLを侵食時間（遅れ時間）で割ったものとして捉えられています。遅れ時間は粒径dmが大きければ長く、流速 U が速ければ短くなる形として表現されており、補正係数αで調整されています。この表現が正しいかどうかはわかりませんが、αを合わせこみパラメータとしていますので、ある程度の調整はできるのでしょう。

2013年 Hyper KANAKO 勉強会の第一部で紹介されたスライドによれば、この式をさらに簡略化して実装しているようで、その影響をαに含めています。

堆積速度も遅れ時間を考慮し、βといった合わせこみパラメータを使っています。

つまり、支配方程式の一つ、最も重要な河床変動の式 ∂z/∂t=i が「高橋のモデル」では合わせこみパラメータとみなせる扱いです。先日の河床変動計算では理論式でしたが、これは驚愕です（権威の届かない海外では通用しないでしょう）。

なお、最大粒径10㎜での室内実験では、α=0.0007（遅れ時間5秒程度）が示されています。またそれを参考に、再現計算ではα=β=0.0001を用いています。で、最終的にはβとαの一般性についての検討は今後の課題とされています。

一方、Hyper KANAKO のスライドでは、式が文献に比べ簡略化にもかかわらずα=0.0007が良く使われるとコメントされています。マニュアルでは以下の通り。大丈夫なのでしょうか。ま、結果が異なるのは、合わせこみパラメータを調整せずにデフォルト値を使うユーザーが悪いとも言えますが。

Coefficient of erosion velocity 侵食速度係数 δe 0.0007

Coefficient of deposition velocity 堆積速度係数δd 0.05

δe・δd は、高橋モデルで日本の石礫型土石流を計算する場合は、この表に示した値を用いる例が多い。ただし、上手く合わない場合は、これらをチューニングする方法が、よく取られる。

さらに文献を探してみると、式の補足的説明が以下の講座に掲載されていました。※式22：侵食速度式、式23：堆積速度式

江頭「土石流の停止・堆積のメカニズム(1)」砂防学会誌46巻(1993)1号

https://www.jstage.jst.go.jp/article/sabo1973/46/1/46_45/_article/-char/ja

式(22),(23)あるいは式(23)は見かけはすっきりしているようにみえますが,これらを直接検証するようなことが行れていないことも原因のひとつですし,最近, 共同で行った実験(寺西・江頭,1993)によれば,堆積速度は,式(22)の傾向とは逆にh/dが小さくなるに伴い増加するような結果が得られている

27年前でこの知見です。それから現在まで何も変わっていないのは、検証済みなのか、それとも放置なのか。
いずれにしても、再現計算なしでは非常に危ういモデルであることは理解できました。

土石流モデル 不明点の整理

久しぶりに Hyper KANAKO を使う必要が出てきたため、この週末にわかっていること、わかっていないことの整理を始めました。

このソフト、理論マニュアルがないため、どのような式が使われているのか一見わかりません。マニュアルがないため文献を追いかけることもできません。開発者に問い合わせることで論文等を紹介されたことがありましたが、そこにも詳細は書かれておらず、最終的には引用のループでお手上げ状態でした。

ソースが公開されていないため、そちらから追いかけることもできません（後輩さん曰く、'嫌気が差すほど’不具合が多いのですが、残念ながらその改良もできません）。

https://phreeqc.blogspot.com/search/label/KANAKO

このような状況なのですが、土研さんの資料4240号がこのソフトに準じておりにソフトが適合しており、特許や著作権のしがらみがなく使用できる唯一のソフトであるため、研究発表会などでは時折耳にします。開発が学主体なのでお客様にとっては「お墨付き」があるように受け止められるのも一因でしょう。

で、先輩さん曰く「これを使わざるを得ない」と。

で、現段階での整理です。

わかっていること

• 支配方程式4種
  1.運動方程式、2.３.連続の式（水＋土砂、土砂のみ）https://phreeqc.blogspot.com/2014/12/blog-post_6.html
  ４.河床変動の式　∂z/∂t=i 

• 侵食・堆積の判断、およびその速度 i の算出に平衡土砂濃度（土石流の段波が一定高さを保つ流下土砂濃度=流下中に侵食も堆積も発生しない濃度）を利用。
• 平衡土砂濃度は、流下土塊＋河床土塊を一体（剛体）とみなした、ある深度でのせん断応力とせん断抵抗力の比より導出（高橋「土石流の発生と流動に関する研究」）
• 平衡土砂濃度は水面勾配（エネルギー勾配や河床勾配より安定のため代用）を用いて算出。
• 3つの形態（土石流、掃流状集合流動、掃流砂）に区分し、それぞれ平衡土砂濃度の計算式を変更。シームレスに切り替え。
• 河床せん断力τの計算式も、河床土砂濃度と流下土砂濃度に応じて3区分

わかっていないこと

• 河床せん断力τを３区分する濃度の閾値（0.01、0.4C*）の根拠
• 河床せん断力τの導出（土石流、掃流状集合流動）
• 侵食堆積速度 i の導出（特に侵食・堆積速度係数の一般値の根拠）
• 平衡土砂濃度c∞を３区分する水面勾配の閾値（0.03、0.138）の根拠
• 平衡土砂濃度c∞の導出（掃流状集合流動、掃流砂）

要は土石流特有の「高橋らのモデル？」※ を追えないということです。

実験で決めました、としか言いようがないのでしょうか？

*********************************

※「江頭モデル」との比較をおこなっています。

高濱ほか2000「土石流シミュレーショモデルの比較研究」https://www.jstage.jst.go.jp/article/prohe1990/44/0/44_0_687/_pdf 

*********************************

20200602追記。

2019年版の Hyper KANAKO も不具合が多く、1日頑張りましたが計算までたどり着けませんでした。実務には向かないでしょう。

また寝かせましょう。

Hyper KANAKO と Civl3D

年末のテコ入れで、それなりに動く様になった Hyper KANAKO。

後輩が頑張って計算をかけているのですが、まだ不具合が出ている様です。

計算・書き出し・ファイル容量を軽くするため、一次元領域を 50 mピッチに変更していたのですが、それでは結果の断面表示に問題の生じることがわかりました。ま、計算値はテキストで書き出されていますから、なんとかなるでしょう。

天然ダムの縦断方向の延長は入力できませんので、一次元のピッチに合わせて複数の入力になります。案外、その方が利点もありました。天然ダムの横断幅と高さを個々に反映できるのです。

KANAKO での天然ダム位置入力は GIS 上でのマウス指定ですので、複数の横断箇所を一定ピッチで入力するのは困難です。Civil3D で縦横断を切り、帯の値などを EXCEL に手入力したほうが容易です。

幸い、KANAKO では流路が SHP ファイルになっていますので、それを Civil3D で読めば、線形に指定できます。（と思ったのですが、あまり綺麗な形状で読めなかったので、ArcGIS 経由で Polyline 化し、読みました。）

崩壊前地形、崩壊後地形、線形を利用し、Civil3D で縦横断を作成。それを後輩に返し、KANAKO に反映させ、計算！

うまくいったようです。

KANAKO の次 Ver. 公開が、今月上旬というアナウンスでした。制限や不具合が取れていることに期待しましょう。

Hyper KANAKOの制限

先週より、「KANAKOが動かない」と後輩に言われておりました。

先輩のお手伝いが終わり、今度は後輩の番。私もHyper KANAKO に詳しいわけではないのですが、待たせた分、答えを出さねばなりません。できれば短時間で。

まずは何が動かないかの検証です。
後輩の作ったモデルはハイドロが数日、1次元が数km、天然ダム多数、といった予想よりかなり大きなものでした。おそらくモデルの大きさ、確保したメモリ−関連で「動かない」と言うのでしょう（このような時、ソースを確認できない点はネックです）。

問題を切り分けるため、LPデータはそのまま、1次元、2次元範囲は縮小、ハイドロはデフォルト、定数や計算条件もデフォルトに戻して計算。結果、うまく流れます。

次に天然ダムを1つ入れて計算。これも決壊し、流下します。

次にハイドロを数日、計算時間を 50000秒に設定。で、ここで OUT。計算途中でエラーが出て最後まで計算しません。ハイドロの土砂濃度を低くし、流れやすい設定に変えてもダメ。同じ時間で止まります。

刻みを 0.01秒から 0.1秒に変更すればその 10倍の時間まで計算します。これは1次元・2次元の計算領域、回数、確保したメモリー、書きだすファイルの大きさなどの制約だろうと踏んで、領域を変えてみることに。

結果、2次元領域を変えても変化なし。
一方、1次元領域を2.5倍にすると、計算が止まる時間が約1/2.5に減りました。1次元の計算は5mピッチ固定なのですが、これを手入力で 50m にしてやると、計算が止まるまでの時間は 10倍に増えました。ファイル書き出しを含む1次元の計算に制約があるようです。

対応としては、開発元である大学の先生に問題を報告・確認し、1次元の制約であれば改良を申し入れるというのがBESTです。が、すでに冬休みでしょう。
対処療法ですが、1次元の範囲を可能な限り狭く、1次元の仮定を逸脱しないピッチと刻みを使って制限内に収めて乗り切る、といったところでしょうか。「制限」は、そのモデルでのアタリ計算からある程度導けるでしょう。

ま、これで先輩、後輩のお手伝いは終わり。答えも出せて一安心。
次は、自分の番です。

問題解決の視点

天然ダムの決壊計算を実施し、「まだまだ」と思った点に以下のようなものがありました。

・浸透流やパイピングの取り扱いがない。

あくまで手を動かした KANAKO での感想ですが、理論上、他の多くの河床変動計算コードにも当てはまると思います。

天然ダムを越流するまでは、計算上、ダムが安定しています。そこに浸透流やパイピングの概念は組み込まれていません。流水による河床変動の計算をベースとしているので当たり前なのですが、実際とは異なり違和感満載ですね。

話は変わりますが、昨年発生した広島の土石流について、地下水の上昇のみでは説明できない点がいくつか挙げられています。

・1次元浸透流で雨を入れると、表層部は飽和するが底部は飽和しない。

→降雨が短時間のため、空気が封入され飽和しない。

・室内試験による強度定数を使用すると、安全率が1.2を上回る。

→ではなぜ崩壊したのか？

・軟岩が削られ、V字・Ｕ字の谷が形成されている。

→6.29ではなかった現象。

被圧地下水、パイプ流、空気圧が関与した可能性があるようです。深層崩壊と共する部分があります。それらを計算に載せるためには、綿密な踏査による地質状況の把握と調査法・解析法の技術開発といった両輪・フィードバックが必要でしょう。

問題を解くには既存の幾つかの視点、新たな視点の組み合わせが重要になるということだと思います。一歩引けば、見えてくるものもあるのでしょうね。

*************************************
20150811追記
同僚や先輩が、昨年度発生した広島土石流に関わっています。
相談を受けてアドバイスしていましたし、講演も聞いていましたので知ったつもりになっていました。報道も可視化結果も見ましたし、6.29にも携わっていましたので、余計に勘違いしていたのでしょう。

先日、現場を眺める機会があり、そこで勘違いに気づかされました。

想像より大きな現場でした。机上では伝わらないものがあります。

今後に活かすためにも、現場と机上、両立しなければと反省。現場を見るだけでもダメ、可視化やシミュレーションなど、机上だけでもダメ。両立できて初めて人の役に立てるのでしょう。

Hyper KANAKO その2

Hyper KANAKO の天然ダム決壊解析機能を試してみました。

最初は2次元領域にダムを設定し、エラーの続出。
そういえば、説明資料に1次元計算であることが書かれていたと思い出し、1次元領域に配置し直して再計算。これで最後まで流れました。

で、結果の表示になりストップ。
描きたい断面の位置を手書きで指定し、「結果断面図作画」ボタンを押してみるのですが、無反応。
1時間ほど考え、ようやく原因が判明。どうやら、計算結果は MapWindow のプロジェクトフォルダでなく、デフォルトでは以下に保存されるようでした。これを指定しないと、読み出せなかったわけです。
C:\HyperKANAKO\ThreeD\Simulation\（プロジェクト名）\

設定を正しますと、断面の表示が確認できました（時々エラーが発生しますが）。
高さ40mの天然ダムを設定したのですが、どうも河川縦断方向の延長は指定できないようです。従って極端に薄いダムができてしまいました。天然ダムが薄いため、計算結果も決壊以降の下流側への土砂供給が少なすぎ。実際とは全く異なる形になりました。

操作は簡単ですので、もう少し開発が進み、安定性が増せば、実務でも利用できそうです。それまで再度、眠らせておきましょう。

Hyper KANAKO

講演に影響を受けたこともあり、ずっと眠らせていた Hyper KANAKO をひっぱり出して、手を動かしてみました。

気にはなりつつ、実験式を追えない（出典に詳細が書かれていない）ので、机の中で眠らせていました。

HPから最新版をインストールしてUSBを更新しセットアップ。一通りのマニュアルを読み、サンプルファイルを動かす。これだけなのですが、案外あっさり動いて、昼までに終了。堰堤の入力や計算結果の表示も問題なく動きました。

昼から広島の土石流発生箇所を題材として計算してみました。
基盤地図情報の10mメッシュを変換するツールがあったので、5mでなく10mをDL。それを変換後、初期設定ファイル作成まではあっさり進みました。
が、1次元の計算範囲を決める際にストップ。GUI は MapWindow ベースなのですが、Web 上の地図が取り込めません。OpenStreetMap 等も読めるように開発側は考えているようですが、完全対応はVer.4.9からのようです。Web Map Service (WMS) 等が容易に利用できるようになっていますので、数値地図25000は久しく使用していません。が、当面はそのCDが必要なようです。レトロですね。残っているか探してみましょう。

また、計算結果では、なぜかゴミが発生。1次元の計算部分から2次元平面領域に入った直後、計算領域の端の方で水位が発生しました。何度か実施してみましたが、すべてに発生します。どこか、設定を間違えているのでしょう。

VRMLの表示にも問題。計算結果を段彩図として重ね、VRML で 3D 表示できる仕様のようですが、それがうまく動作しません。GIS 内では表示できますので計算は（ゴミが出ていますが）正常に終了していると思われます。VRML 用のファイルを作成する段階で、何かエラーが発生しているようです。

何度か試してみましたが結果は同じ。結局、解決しませんでした。

ま、急ぐ話ではないので、次に天然ダムを取り扱ってから、また寝かせておきましょう。

連続の式・運動方程式

まずは大雑把ですが、土石流（KANAKO）の「連続の式」と「運動方程式」をやっつけましょう。

参考にした資料は以下の通りです。

• 河村哲也「河川の流れのシミュレーション」山海堂 pp.1-23（導出はこれ）
• GUIを実装した土石流一次元シミュレータ開発：中谷加奈，里深好文，水山高久，砂防学会誌，Vol.60, No.2，pp.41-46, 2008
• 土石流計算における1次元・2次元シミュレーションモデルの結合，和田孝志，里深好文，水山高久，砂防学会誌，Vol.61，No.2， p.36-40，2008
• 中村ほか「地震砂防」pp.121-123 
• 土研資料「深層崩壊に起因する土石流の流下・氾濫計算マニュアル」p23（KANAKO 使用前提のマニュでしょうか？それとも砂防では基本的な式なのでしょうか？）

＊以降、誤りがあれば、随時修正します

～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

ρを定数とし、鉛直方向に積分した連続の式は、以下の形になります。

∂h/∂t＋∂M/∂x = a

ただし、M = Uh　（流束）で、U は u の深度方向の平均速度です。

LSFLOW では a の部分が 0ですが、KANAKO では河床変化につりあう水＋土砂量（河床変動の式侵食・堆積速度）となっています。また、別途、土砂濃度を両辺に乗じ土砂量のみのつり合いも別式として加えています。
ちなみに浸透流では鉛直方向に積分する必要がないので M が vi になり、ρの変化がないと仮定すれば第1項が消え、右辺のaが注水・揚水、地盤が保持する量として考慮されます。

～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

運動方程式は、LSFLOW と KANAKO が同類ですね（浸透流はダルシー則）。
どちらもナビエ・ストークスを鉛直方向に積分して得られる「河川の基礎方程式」のようです。ただし、KANAKO は初めから粘性項を無視しています（下記の通り）。LSFLOW では、実務上粘性を0とみなす事が多いので、実質同じでしょう。uが深度方向に一定 = U と仮定すると、以下の通りとなります。

∂M/∂t＋∂uM/∂x= -gh∂H/∂x-τ/ρ

この中の∂H/∂xを水面勾配、τを河床でのせん断力で表現しているのが KANAKO。ただし、せん断力の中身については、まだ理解できていませんので、ごっそり後回しです。
LSFLOW はτをクーロン則、ニュートン型の粘性抵抗則、マニングの抵抗則として選択できます。

～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

20200530一部修正

土石流の支配方程式

Hyper KANAKO をさわることになりました。

http://hyper-kanako.com/

昨年の LSFLOW が比較的容易でしたので、今回も、同じ平面2次元であれば、何とかなるだろうと踏んでいます。LSFLOW は既に岩盤のほぐし率まで考慮できるよう改変済みなのですが、今回、私の頭が追い付けば KANAKO も統合してみたいと目論んでいます。

が、甘くないですね。早々躓きました。土石流の支配方程式が理解できません。

運動方程式は比較的スタンダードな形ですが、x方向のせん断力をh（水深）で割っている理由がイマイチ分かりません。hを定数とみなし、全体をhで割った結果でしょうか？連続の式も右辺が0でなく、数値が入るようです。ま、これは河床変動とつりあわせるための修正でしょう。

一番とっつきにくいのが濃度。土石流の土砂移動は濃度で扱うようです。これには驚きました。

平行濃度とは何ぞや？土石流・掃流状集合流動・掃流砂の濃度閾値の意味は？その式の中身は、など疑問点山積みです。ただ、先輩に聞けば、河川砂防の分野では常識的だそうです。確かに、どの論文をみても導出がないので常識的なのでしょう。運動方程式や連続の式はすぐ追えると思いますが、濃度は自信ないですね。

まだ時間はありますので、一つずつクリアーして行きましょう。