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

まずは、侵食/堆積速度 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日頑張りましたが計算までたどり着けませんでした。実務には向かないでしょう。

また寝かせましょう。

GSFLOW-GRASS と ModelBuilder

2年前はココまででした。

GSFLOW-GRASS（途中）

https://phreeqc.blogspot.com/2018/08/gsflow-grass.html

今回は Win10 に GRASS GIS を入れ、

install_extensions して、

python domain_builder\buildDomainGRASS.py DRIVE_LETTER:\PATH\TO\settings.ini するとエラー。

GRASS が Python を認識しません。うーん。

Python 周りの環境を触りたくないので、Ubuntuへ。

Ubuntu なら Docker でしょう、ということで、docker-compose build。

で、エラー。

2年前の構成なので、リンク先のアーカイブ等が更新・変更されていました。

それらを書き換えてエラーを修復し、再度 build。

で、docker-compose up、docker exec -it CONTAINER ID bash、grass。

が、grass を認識しません。うーん、これ以上は実力不足。

以前作っていたイメージを消してしまいましたし、野良イメージの中身を確認するのも面倒なので、Ubuntu もあきらめました。

で、基本に戻り Win10＋ArcGIS Pro 2.5。vegetation 情報を与える直前まで実施。

手作業だと、着実にデータを作成できるのですが、遅い。ま、時の流れと共に、ArcHydro や XTools が不要になっているのはありがたいことです。

ふと、ジオプロセッシングの履歴が画面右側に蓄積されていくことに気付きました。

これ、保存できないでしょうか？

サポートさんに伺うと、プロジェクトファイル毎に履歴が保存され、しかも Model Builder に D&D できるとのこと。

これはすばらしい！と早速やってみると、できました。もう一度、素晴らしい！

といっても、今回は途中までの作業でしたし、所々に重複した履歴が入っていましたのでクリーニングが必要でした。が、これがあれば当面、GSFLOW-GRASS を稼動できなくても支障ありません。

GSFLOW-GRASS に関しては以前より後退。

ArcGIS Pro では大幅前進。

次はどうなるでしょう。

河川砂防技術基準の河床変動計算

HEC-RAS の1次元河床変動計算が、河川砂防技術基準に適合するかの確認です。

調査編 第6章に該当するのですが、恥ずかしながら詳細を追ったことがなかったので、この機会に頭に入れました。

河川砂防技術基準　調査編（平成 26 年 4 月）

第 6 章 河床変動、河床材料変化及び土砂流送の解析

https://www.mlit.go.jp/river/shishin_guideline/gijutsu/gijutsukijunn/chousa/pdf/05.pdf

1．1 位置付けと目的

沖積河川を対象

沖積平野を形成する河川という意味でしょう。直轄区間を想定されているのでしょうね。そのため、山地河川は対象外と読めます。準じて使用といったところでしょうか。

1．3 解析レベル、第2節 目的に応じた解析レベルの設定

１）河床高の解析レベル：1 次元[1DB]、2 次元[2DB]

　　平均河床高：[1DB]

　　局所洗掘深の予測、河床高の平面分布：[2DB]

２）粒径の解析レベル：一様粒径[U]、混合粒径[M]

　　対象区間の全域で材料の粒径範囲が同一→一様性が高い

　　粒度分布の淘汰がよく、粒径範囲が狭い範囲に集中している→均一性が高い

　　淘汰が悪く、広い粒径範囲にわたって偏りない粒度分布→均一性が低い

　　一様性と均一性が共に高い場合：[U]

　　それ以外：[M]

３）流砂の解析レベル：掃流砂[BL]、掃流砂と浮遊砂[BSL]

　　礫：[BL]

　　砂又は 砂礫：[BSL]

４）流れ場の解析レベル：1次元[1DF]、準2次元[2DF’]、2次元[2DF]、2次元（2次流付加）[2DF+]、準3次元[3DF’]、3次元[3DF]

　　河床高の解析レベルと対応

　　[1DB]：単断面[1DF] と複断面[2DF’]

　　[2DB]：[2DF]～[3DF]

1．4 解析の実施手順 

１）解析レベルの設定

２）解析法の設定：計算法と各種条件等の設定

３）解析の実施と検証

第3節 計算法の設定 

3．2 流れの計算法の設定 ＜標 準＞ 

[1DF][2DF’][2DF][3DF’]については、第5章河川における洪水流の水理解析の 4.1、4.2、4.3、4.4、4.5、4.6に基づいて設定することを標準とする。

第5章4.1は縦断におけるエネルギー（energy grade line）の保存式です。4.３以降はこれをベースにしています。先日組んだ「現場のための水理学」やHEC-RASも同じ考え方です。単純なのでこの分野では共通なのでしょう（厳密には、基準では運動エネルギーの係数1/2が示されていません。導出が異なるのでしょうか？）。

3．3 河床形状の計算法の設定 

河床形状ｚBは流砂の連続方程式により算定する。 種別【1B】【2B】【1BS】【2BS】の流砂の連続方程式は式(6-3-1)～(6-3-4)を標準とする。

掃流式(6-3-1)は「現場のための水理学」式5.3と同一。HEC-RASも同一です。6-3-3は6-3-1に浮遊砂を足しています（後述の代替え案を使うことで、結局6-3-1式と等価になります）。混合粒径に適用する場合のパラメータは3.5~3.7で示す計算法によるとのこと。

3．4 粒度分布の計算法の設定 

交換層の枠組みを用いるのを標準。

粒径範囲が広い河床材料、2 峰性の粒度分布を有する河床材料、及び石礫で構成される河床材料 等への適用性（各種の流砂量式を含む）も適宜参考として検討が重要。

河床の各粒径階の含有率 fbi は、各粒径階の連続方程式により算定する。種別【1fB】【2fB】【1fBS】【2fBS】のfi連続方程式は式(6-3-5)～(6-3-8)を標準とする。

これも混合粒径、浮遊砂に適用する場合は3.5.2、3.7で示す計算法によるとのこと。

HEC-RAS の hydraulic reference manual では、ここで初めて日本人の名前を目にしました。基準でも基本的考え方の資料1)として紹介されている文献です。先駆的研究だったのでしょう。文献内で式(6-3-5)が示されています。

https://www.jstage.jst.go.jp/article/jscej1969/1971/195/1971_195_55/_pdf

この手法は、最もシンプルな Active layer Mixing Method として実装されています（アーマリング因子がこの手法には含まれていないため、利用に注意とのこと）。

3．5 掃流砂の計算法の設定

3．5．1 一様粒径の掃流砂 （１） 一次元解析の掃流砂量式 

例示として[BABL]の掃流砂式が挙げられています。式の設定は「決め」ではなく、検証・調整項目の一つと位置付けられています。

3．5．2 混合粒径の掃流砂

交換層 における粒径階 diの含有率 fbiを乗じた式(6-3-21)を用いるのを標準とする。

 qBi=fBif(di,τ,τci,ρ,σ,g,etc )

考え方として示されており、○○式を使用といったような規定はありません。

hydraulic reference manual 12-47式、12-54式でも同様の表現が使われています。

3．7 土砂濃度分布の計算法の設定

土砂濃度分布 C は浮遊土砂の移流拡散方程式により算定する。 

これは驚きました。今まで見たソースでは、このような計算はなされていません。と思ったら、次の例示中に代替え案が用意されていました。

河床縦断形状の長期変化等の土砂濃度分布の非平衡性が無視できる事 象を対象とした場合には、土砂濃度分布 C を式(6-3-35)により計算するのに換えて式(6-3-42)により算定し、これを式（6-3-3）、(6-3-7)に代入して河床高及び粒度分布を計算す ることができる。 

なるほど。濃度は使わなくてもOKです。ほとんどのソフトが移流拡散を解いていないのではないでしょうか？

結局、これを使うことで「3．3 河床形状の計算法の設定 」の式6-3-1と同じ形になります。

HEC-RASの場合、「3.5.2」の段階で掃流砂 or 掃流砂量＋浮遊砂量を扱える式を選択できるようになっています（量でなく、すべて濃度として扱う式もあります）。「3．3 河床形状の計算法の設定 」の段階で両者を合わせるのではなく、その代入前に合わせた量を求めておくイメージでしょう。

第5節 解析の実施と検証・調整 

5．1 解析の実施 

流れ、流砂、河床状況の計算は、時間項について逐次数値積分することによって行い、解を得ることを標準とする。 その際、流れに対して河床の時間的変化は一般に遅いため、流れを解く際には各時間ステッ プ内での河床の時間的変化は無視する「準定常」の条件で計算を行う。 

Quasi-unstedy のことです。疑似非定常と訳していましたが、これからは「準定常」を使いましょうか（「準定常」は河川屋さん以外には通じないでしょう）。

5．2 解析の検証・調整

調整は、粗度係数の設定、交換層厚の設定（3.4）、各種流砂量式の設定(3.5、3.6、3.7)、及 び初期条件・境界条件等の設定(第 4 節)が対象となる。 流砂式の補正（パラメータの補正等）も必要に応じて調整に含めてよい。

これは自由度が高いですね。

以上、再現結果で精度を確保できるなら、ある程度自由に計算式等を選択しても良いといった内容でした。性能重視、技術者の力量を問う、といったところでしょうか。

HEC-RAS は問題なく基準に適合していることが分かりました。日本人の提案式が採用されていませんのでお客様の反応が気になるところですが、基準上は問題ありません（前に聞いたガラパゴスとはこのことだったのでしょうか？）。

国内で流行っていない理由は、技術的内容とは別のところにあるのでしょう。

HEC-RAS Ver.5.0.7

HEC-RAS には、以下の機能がありました。

・1次元定常で常流、射流、混在を選択可能。

・1次元非定常、疑似非定常で河床変動計算可能。

・掃流砂、浮遊砂を扱える（しかし、日本人考案の式は出てこない）。

・粒度を扱える。

・支川上流端で非定常、疑似非定常の土砂流入条件を設定できる。

・複数の支川の合流、分岐を扱える。

十分です。

基準との適合性を確認してみましょう。

HEC-RAS vs Mike11

Fully conservative coupling of HEC-RAS with MODFLOW to simulate stream–aquifer interactions in a drainage basin 

https://www.sciencedirect.com/science/article/pii/S0022169408000772

どうやって動かしているのかわかりませんでした。非定常で計算できる優位性があるそうですがイマイチしっくり来ていません。

先日、HEC-RAS を触っていたら、メニュー内に MODFLOW の文字がありました。グレー-アウトして使えず、マニュアルにも含まれていません。開発中の機能が文献に出てきたのかもしれません。

Comparison between capabilities of HEC-RAS and MIKE11 hydraulic models in river flood risk modeling (a case study of Sungai Kayu Ara River basin, Malaysia)

https://www.inderscienceonline.com/doi/abs/10.1504/IJHST.2012.049187

こちらは MIKE11 よりも精度が良かったとあります。なら、MODFLOW の相手として選択すべきでしょうか？ MIKE11 + MODFLOW は市販されていますから、ユーザビリティの点では Mike11 の方が一歩リードですかね。

これらの組み合わせ、GSFLOW より使い勝手が悪そうですが、モデル化は容易でしょう。HEC-RAS 自体、高機能であることが理解できつつありますので、手放さない可能性大です。Ver.UPの度に連成機能をチェックすることにしましょう。

Servant or Master

Application of classification trees for improving optical identification of common opaque minerals

https://www.sciencedirect.com/science/article/pii/S0098300418310604

決定木を用いた鏡下での不透明鉱物の分類です。

決定木がどのように判断・分離したかを、人が解釈できます。が、完璧に区分しようとして多岐に渡ると煩雑になります。

この文献では、人が判断しやすいように3つの特徴量のみで答えを出すように加工されています。機械学習の王道とも思える利用法です。

機械学習で多種多様、多量のデータを人が解釈しやすいようにうまく加工する。その結果を用いて、これまで通り技術者が何らかの意思決定を行う。扱える情報量が格段に増えた現状において必要とされる裏方技術が増えただけです（今月の実践ビジネス英語では「servants working in the background」と表現されています）。

このような裏方作業の結果を、さらに機械が判断するようにもなりつつあります。製品化目前の自動運転などがその例です。人と共存する技術に高まりつつあります。

もちろん、人がスキルを放棄しないことも重要です。機械とは一線を引き、利用するが依存しない状態でとどまることも必用でしょう。

何かしら自分を助けてくれる Servant を手元に置いておくことは今後必須です。その時、利用者としての Master の能力がより問われるのでしょう。

将来、どちらが Master になれるかの争いは、既に始まっています。

地域特化と過学習

扱う災害発生降雨が1降雨のみであれば、非発生降雨と完全に分離できます。

数十件程度までなら区分できるでしょう。単なる過学習ですが。

昨日書きとめた文献のように、自治体単位だと過去の発生事例が限られるため、一見、識別性能の良く見えるモデルができるでしょう。が、当たるのか？と言われると、おそらくダメ。同条件で同じ降雨が来ない限り、このモデルでの予測が当たる保証はありません。機械学習分野では避けなければならない考え方です。

では、全国レベルで学習して地域が当たるのか？と言われると、個人的にはある程度当たると考えています。経産省のようにデータを公開してコンペに出せば、すぐにデータサイエンティストたちが一定のレベルまで引き上げてくれるでしょう。解決策の一つです。

地域のデータを深堀りすることは必要です。が、将来に役立てるためには、より大きな視点に基づいた評価が必要です。地域特化が多量のデータを扱える技術がないことの隠れ蓑になってはいけません。過学習に陥ってはいけないのです。

全国かつ非発生にも目を向け、恣意的にならないよう、繰り返し検討を進める必要があります。

再現期間と崩壊

貯めていた「再現期間と崩壊」に関する文献をチェックしました。

齋藤ほか「日本列島における斜面崩壊発生と確率雨量との関係」

斜面崩壊は、その場所の暖候期降水量にかかわらず、再現期間 5 年（中央値）～ 300 年（95th percentile）の降雨イベントで発生していた．

その中でも規模が 1千 m3 以下の斜面崩壊は，再現期間 10 年以下の降雨イベントで発生。

10万 m3 以上の斜面崩壊が発生した降雨の再現期間は 100 年以上であった。

木下ほか「深層崩壊が集中的に発生する降雨条件―平成２３年台風１２号の降雨分析―」

台風１２号の降雨は短期雨量指標ほど超過確率年が小さく、長期雨量指標ほど超過確率年が大きい降雨イベントであった。

各降雨指標と深層崩壊発生の関係には明瞭な相関関係が見られなかったが、48時間雨量および累積雨量の確率年と深層崩壊の発生数および発生率、深層崩壊規模の間には相関関係が見られる。

 特に崩壊面積10万m２以上の崩壊については累積雨量１/180年以上、48時間雨量１/120年以上で発生する傾向が見られる。

広島の土石流が約10年毎の発生ですので、再現期間10年というのは妥当な気がします。

また、深層崩壊など大規模なモノで100年以上というのも感覚的に受け入れられます。

執印ほか「確率年を用いた単一降雨指標による土砂災害発生危険度の評価について-東京都伊豆大島を事例として-」

IRPIによる土砂災害発生時の下限値は104年、非発生豪雨の上限は64年。

IRPI:スネークライン原点からの距離（確率年）

IRPI は忘れていたのですが、私も2軸のみならず3軸まで試したことがあります。が、全国規模では1軸に比べ特別に効く指標とはなり得ませんでした。

この文献は他の文献とは異なり、非発生降雨に触れられているところがgood!

どのような問題でも、偏った視点だけでは解けません。 

川越ほか「降雨極値の再現期間を用いた斜面災害の発生確率」

1991年9月19日のいわき市での降雨イベントを利用。浸透計算（動水勾配）、起伏量を特徴量としたロジスティック回帰による崩壊確率モデルを構築。

再現期間30年の降雨極値により東北地方の概ねの山地、丘陵が発生確率80%以上になる。

再現期間5年の降雨極値により釜石市、いわき市、再現期間30年の降雨極値により秋田市、仙台市、郡山市が発生確率80%以上になる。

川越ほか「温暖化に対する土砂災害の影響評価 」

2004年7月12日～13日の新潟・福島豪雨イベントでモデルを構築。

再現期間30年などで、山頂から山裾側にリスクの高い地域が拡大することが明らかにされた。※発生確率主体の部分省略

多重ロジスティック回帰は、単純な割に案外良い予測性能を出す印象があります(GBDTより劣るのでモデル採用には至りませんが)。
残念ながらイベント数が限られており、モデルの予測性能も不明です。

いずれにしても、全国の降雨を同一指標で扱おうとして再現期間、確率年を利用する方法に至るのは当然の流れでしょう（全国で指標を統一すべきかはわかりませんが）。

RBFN値と共に再現期間等の利用が進むのか、あるいはビッグデータを扱う新たな流れに呑まれるのか。呑まれたとしても、このような地道な評価、解釈は続けて頂きたいものです。

Dam Break using HEC-RAS

HEC-RAS を利用した Dam Break の文献です。

Landslide dam and subsequent dam-break flood estimation using HEC-RAS model in Northern Pakistan」

https://link.springer.com/article/10.1007/s11069-012-0361-8

A Dam Break Analysis Using HEC-RAS

https://ascelibrary.org/doi/10.1061/9780784482353.036

いずれも、河床変動は扱っておらず、水位のみ検討対象とされています。破堤時の洪水といった問題ですので、それでよいのでしょう。それよりも、HEC-RAS で Dam Break を扱える点に驚きました。

HEC-RAS を利用した文献は、海外で数多く引っ掛かりました（国内では引っ掛かりません。砂防屋さん曰く、日本はガラパゴス化しているとのこと）。支川からの流入や境界条件としてunsteady の土砂流入も扱えますので、利用されている方は多いのでしょう。Civil3D だけでなく、ArcGIS のアドオン？もありますから、環境は整っています。

日本でも CAD や GIS に組み込まれるほどの良いソフトが出てきてほしいものです。

River and Flood Analysis 2020

River and Flood Analysis が動くようになったのかを確認。

2018ではダメでしたが、2020ではどうでしょう？ということでインストール。
https://phreeqc.blogspot.com/2018/05/river-and-flood-analysis-2018.html

無事、メニューも出て（なぜか2回目は隠れましたが）、ボタンを押せば起動します。

早速ためそうと思いきや、これ、河川として Civil3D の線形を指定できません。ポリラインのみ指定可能です。とすると、図面とモデルで測点が微妙にずれるかもしれません。
そうであれば、Civil3D 標準の[出力]タブ  [書き出し]パネル  [HEC RAS に書き出し] からジオメトリを書き出して、HEC-RAS に持って行った方がリスクを負わず効率的です。
https://knowledge.autodesk.com/ja/support/civil-3d/learn-explore/caas/CloudHelp/cloudhelp/2020/JPN/Civil3D-UserGuide/files/GUID-6E2DFC49-041A-4661-8A97-7776A7C16A8F-htm.html
つまり、線形・測点がある場合は標準機能、それらのない場合は River and Flood Analysis モジュールを選択すべきということでしょう。

ポリラインを作って縦断として認識させましたが、ココでも失敗。

下流からポリラインを書き始めたのですが、Ｒ&F で認識させると書き終わった位置が起点になります。上流から書き始めるか、事前にポリライン編集で反転しておかなくてはなりません。ま、このあたりは慣れでしょう。

縦断を認識させると自動で横断が切られ、右岸左岸をポリラインで指定しておけば全横断通してマニングの粗度係数を適用でき、流量を設定したらあとは計算だけ。これはかなり効率的です。横断をテキストデータにしなくても計算できてしまうのはありがたいですね。

最後に結果の描画。水面を地形に合わせて描いてくれます。Ｒ&F ではこの機能（水面の自動作画）に期待していました。
が、残念ながら雑すぎてダメでした。作画だけの計算に40分もかかったのに、水位コンターと地形コンターが合っていない、なぜこうなった？というようなところが多々あります。ま、これは横断位置に高さを与えて面を作って切り出しという一般的な作業で綺麗になるので回避できるのですが、現状維持です。

計算用のモデルを HEC-RAS 用に書き出し、不定流（or 疑似不定流）＋河床変動を HEC-RAS で計算させ、結果を Civil3D で作画というのが効率的でしょうか？HEC-RAS の河床変動計算サンプルデータを見ると「研究用」といったようなコメントがあるのですが、どのような扱いなのでしょう？

まだわからないことがありますが、ひとまず計算結果の可視化までは確認できたので良しとしましょう。

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

Civil3D 標準の[HEC RAS に書き出し]では、2020になっても流下方向が逆向きのまま修正されていませんでした。
この問題に対して以前は気づかなかったのですが、以下の方法があったようです。
※メモ書きがありました。単に忘れていただけでした。
https://knowledge.autodesk.com/support/civil-3d/learn-explore/caas/sfdcarticles/sfdcarticles/How-to-reverse-the-stream-direction-in-Hec-Ras-s.html

Issue:
Using "Export to Hec Ras" feature, you export a .GEO file from Civil 3D to Hec Ras. Once you import the .GEO file into Hec Ras, the stream direction is reversed.

Solution:
To Reverse the flow direction in Hec Ras, go to the Geometric Data editor> GIS Tools> Search Invert Lines Table > and Flip Coord. Order.

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

20200605追記

ポリラインを選択時に、以下の表示が出ていたのを見落としていました。

Click on the downstream end of the river centerline:

ポリラインの向きは無関係で、ポリラインの下流側を選択すれば正しく認識してくれました。

Kinematic wave 近似と Diffusion wave 近似

2次元の不等流・不定流計算がどのように導出・実装されているのか調べてみました。

数冊の図書を読んでみましたが、どれも似たような内容でした。そいえばLSFLOW も同じような導出だったなあと思い出しました。LSFLOW のパラメータを工夫すれば地表流の2次元平面の計算ができそうですが、汎用的ではないですね。

そんなことを考えながら新たな図書を読み始めると、惹かれる近似を見かけました。

Diffusion wave （拡散波）近似

関根正人「移動床流れの水理学」p16-19

kinematic wave とあわせて示されています。いずれも扱いやすそうな、より単純な式に変更されています。

忘れていましたが、GSFLOW は kinematic wave 近似を使用しています。地形勾配と逆方向への流れ（バックウォーター）はダメというのが導出を見ると納得です。

最近教えてもらった PRI-model は Diffusion wave 近似でした。地形勾配でなく水位（標高）勾配を使っています。これだと、バックウォーターも計算できそうですが、地形と逆勾配になる地点で引っ掛からないのでしょうか？なぜ'拡散波'と呼ぶのでしょう？知らないことばかりなのですがこの分野では常識なのでしょうか。

いずれの近似も単純なので計算は速いでしょう。タイムステップが全く異なる地下水との連成に使われているというのは納得です。

複断面での疑問点

1次元河床変動計算の続きです。

今日も丸1日時間があったので、単断面から複断面（一般断面）に変更しました。が、よくわからないところが出てきました。

まず、水位の計算方法。昔の河川砂防技術基準には井田法というのがあったそうです。
一方、この手法よりも断面分割法の粗度係数を検討した方が良い結果につながるという文献も見かけました。

30年以上前の「現場のための水理学」では（工夫前の）断面分割法でしたが、現在の主流はどちらなのでしょう？ま、これは河川屋さんに確認すればすんなり解決できると思われます。

次に、断面内での河床を変動させる量と場所。iRIC のマニュアルを見ましたが、詳細は書かれていません。手元のソースでは水面下の河床を一様に上げたり下げたりしています。平均流速に応じ、一様に上げ下げするのはいかにも1次元的です。が、実際は低いところからたまったり、低い部分が集中して削られたりしています。現形状を一様に上げ下げするような変化は起こりません（出水後にまっ平らになった河道も幾度か目にしました）。
計算上、分割ブロックの各流速に応じ河床を上げ下げするという方法も考えられます。ま、これは断差がつきそうですし、そこまでこだわるなら素直に2次元で解けば良いのでしょう。一様な変化とみなしても支障のない問題だけ1次元で解くべきなのでしょう。

他にも、水面下での粗度係数の合成や径深、潤辺の扱い方など、複断面にしただけで多くの疑問点が出てきました。これからも、進むたびに出てくるのでしょう。

ひとまず、ここまで。Quasi unsteady、複断面、複数の支川合流、支川からの土砂供給コントロールまでは、既存コードをトレースすることで実装できました。まだエネルギー勾配の見直しや上記の問題が残っています。一様粒径も簡略化したままなので、きちんと組み直す必要があるでしょう。

この休みの間に手を動かしたおかげで、新たにわかったこと、理解すべき点が明瞭になりました。また、まとまった時間があれば手を付けましょう（忘れそうですが）。

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

20200514追記

HEC-RAS のマニュアルに説明がありました。

水面下の河床を一様に変化させる手法は、Veneer Method と呼ぶようです。やはり、それでは不自然になる場合があり、より深い場所に厚くたまるような実装が加えられているそうです。

メモリ

メモリ192GB のPCで計算。

買った当時、使い切ることはないだろうと考えていたのですが、あっさり使い切りました。最近は 80GB 程度であれば「少なくて済んだ」と思うように。
10年以上前だと 24GB でも奮発したような感覚だったのですが、怖いですね。

GPU だと 40GB で奮発した感覚になっています。10年後だと４00越えでしょうか？

Fortran と Python

Fortran コードを Python に移植し始めて気づいたのですが、Python には goto 文がありません。
調べてみると、可読性の高さを維持するため設けていないそうです。確かに、While 等で代用できますし、あっちこっち飛ばないので追いかけるのが楽です。

地味に面倒だったのが オリジン。Fortran は one-based、Python は zero-based です。+1するだけなのですが、忘れそうになります。

他には ()と[]の表記の違いでしょうか。これは検索・置換である程度修正できるので問題ないでしょう。

Fortran サブルーチンはコンパイルして Pythonで読めるはずですが、今のところ読めていません。これをクリアすれば main の書き換えが主体となります。上記に留意しておけば、ほぼそのまま移植できます。

自動化されないでしょうか？

Python で河床変動計算

GWに入ったものの、出社禁止＆テレワークだと（休日含む）残業を認めないという即席ルールのため、まとまった時間ができてしまいました。

せっかくなので、河床変動計算の復習をすることにしました。簡単な1次元のコードを作りながらより理解を深めようと（あわよくば地表流のコツを掴もうと）いう目的です。

コードは「現場のための水理学」をベースにしました。30年以上前の資料です。今でもBasic や Fortran ソースがHPで配布されています。
http://river.ceri.go.jp/contents/tool/suirigaku.html
まず、これを Python に取り込むことに。メモリ確保の手間や入出力、図化の効率を考えるとこちらの方が圧倒的に楽です。

と、いきなりですが、ここで躓きました。
サブルーチンを dll にしたのですが、Python側で読めません。少し粘りましたが解決せず。時間がもったいないのであきらめて完全移植することに（また後日調べましょう）。速度面のデメリットはありますが、1次元であれば許容できます。

次に改変。
ベースにしたのは不等流であり、これを Quasi-unsteady として時系列データに対応するよう改変しました。条件等をEXCELで作り、pandas 1行で読み込むだけです。本線・支川の流量はもちろん、下流端水位や支川からの土砂流入など、境界条件等を容易に設定できます。ついでにジオメトリもEXCELで。これも読むだけですのでお手軽です。

で、計算。1次元なのですぐに終了。同時に河床形状や水位を時系列で重ねて図化・表示するようにしたので、その場で結果の良し悪しがわかります。この点が特に欲しかった機能でした。

が、その結果がイマイチ。
上流端に境界条件として土砂流入を入れる形にしたのですが、これがマズイようでした。
あらためて上流端付近に支川を合流させ、そこから土砂を入れると良い形状になりました。支川は１つ増えますが、無理のない設定・計算なのでコチラが正解だったのでしょう。1つ賢くなりました。


2・3次元の現象を1次元に落とし込むには計算経験とセンスが必要です。私のような初心者にとって、このような手探りの作業を繰り返すことは必須なのでしょう。
今回、2年前に曖昧にしていた点の一部はクリアになりました。実務で使うにはマダマダですが、今後も改変を続けながら経験を積み上げていきましょう。

TDM-GCC

Fortran を後輩さん（初心者）に使わせたいと相談を受けました。

私は計算速度の点から Intel Parallel Studio を使っているので、フリー版を必要としていません。使うとすれば Cygwin + gfortran や Linux などが選択肢になるのでしょう。

前後の作業効率を考えると、今回は Windows で実装したいところです。しかも手軽に。

探してみるとありました。

TDM-GCC

https://jmeubank.github.io/tdm-gcc/

丁寧な解説はコチラ↓

http://www.compe.media.kyoto-u.ac.jp/syoseki/windows/windows-gfortran.html

手元のソースは問題なくコンパイルできました。計算結果も問題なし。
インストールするだけなので手軽です。Win版でこのような手軽なソフトが配布されているとは知りませんでした。

いくつか試行しないと制限・弱点は見えてこないですが、ひとまずこれが良いでしょう。