LiDARによる土石流観測手法の比較

土研資料とGeomorphology文献の違いをAIさんに整理してもらいました。

基本情報

項目	土木研究所資料	Geomorphology論文
文献名	二次元レーザスキャナによる土石流の流量観測手法	Monitoring debris flow dynamics: insights from 4D-LiDAR observations in Ohya landslide, Central Japan
発行年	2023年11月	2025年
発行元/掲載誌	国立研究開発法人土木研究所資料 第4445号	Geomorphology 482 (2025) 109800
著者	今森直紀、清水武志、池島剛、伊藤誠記	Tatsuki Kaneko, Fumitoshi Imaizumi, Tomoya Osada, Saleh Yousefi, Shoki Takayama
観測地	桜島有村川流域	大谷崩れ(Ohya landslide)、静岡県中部

使用機器

項目	土木研究所資料	Geomorphology論文
レーザスキャナ種類	二次元レーザスキャナ(2D-LiDAR)	4D-LiDAR(3D+時間)
機種	北陽電機 UXM-30LAH-EWA	Livox Horizon (Livox Corp.)
波長	905nm近赤外線レーザ光	記載なし
ステップ角	0.125度	記載なし
測定点数	1520点(190度走査)	記載なし
測距範囲	0.1~30m	記載なし
サンプリングレート	20Hz(50ms毎)	0.1秒間隔(10Hz相当)
測定分解能	1mm	記載なし

観測システム

項目	土木研究所資料	Geomorphology論文
走査次元	1軸走査(断面計測)	3次元走査
設置箇所数	1箇所	2箇所(Site U: 上流、Site D: 下流)
併用機器	非接触型流速計、ビデオカメラ	ビデオカメラ(Sony HDR-CX470、60fps)
電源	UPS経由の商用電源	バッテリー+ソーラーパネル
データ記録装置	ノートPC	Raspberry Pi(オンボードコンピュータ)
起動方式	常時稼働	ワイヤーセンサーによる自動起動
データ保存間隔	連続記録	1時間間隔のファイル保存

計測・解析手法の比較

項目	土木研究所資料	Geomorphology論文
計測対象	流下断面積(水通し断面)	3次元表面形状、縦断・横断形状
座標系	極座標→平面直角座標変換	極座標→直交座標変換(x:鉛直、y:横断、z:流下方向)
解析ソフトウェア	Python 3.9(独自プログラム)	CloudCompare (v2.13)、QGIS (v3.16)
ノイズ除去	1秒間の中央値抽出(雨滴除去)	記載なし
DEM作成	記載なし	複数のグリッドサイズでDEM作成
地形指標	流下断面積のみ	勾配、粗度、これらの平均値・標準偏差

流量算出手法

項目	土木研究所資料	Geomorphology論文
流速計測	非接触型流速計(表面流速×3/5)	ビデオカメラ+LiDAR(直接計測)
流下断面積	LiDARから直接算出(台形公式)	LiDARから3次元的に算出
流量計算式	Qs = U × Ad	詳細な記載なし
時間分解能	1秒(中央値処理後)	0.1秒

観測対象とする土石流の特徴

項目	土木研究所資料	Geomorphology論文
流域面積	約1.6km²	約0.37km²(市ノ沢)
チャンネル特性	砂防堰堤水通し(固定床)	自然渓流(移動床・固定床混在)
河床勾配	記載なし	上流部:2537.3°、下流部:1520°
土石流タイプ分類	記載なし	完全飽和流、部分飽和流
観測イベント数	2022年に4回発生中2回計測	2023年8月3日の1イベント詳細解析

データ処理・出力の比較

項目	土木研究所資料	Geomorphology論文
出力データ形式	CSV(時刻、平均水位、流下断面積、平均流速、流量)	点群データ、DEM、各種地形指標
可視化	ハイドログラフ、断面図、アニメーション(mp4)	縦断・横断プロファイル、地形指標グラフ
プログラム言語	Python 3.9	記載なし(CloudCompareとQGIS使用)
公開データ	プログラムソースコード+計算データ(DVD-R添付)	データ公開の記載なし
ライセンス	プログラムはCC BY-SA 4.0	論文はCC BYライセンス

研究成果・新知見の比較

項目	土木研究所資料	Geomorphology論文
主な成果	・流下断面積の定量的・高時間分解能計測 ・夜間・豪雨時の計測可能性 ・従来手法との比較検証	・完全飽和流と部分飽和流の形態差異 ・表面粗度と乱流の関係 ・堆積過程における逆勾配地形形成
流速に関する知見	表面流速×3/5を平均流速として使用	先端部より後続部の表面流速が高い場合あり
断面形状	台形断面として計算	凸型断面(部分飽和流)の観測
堆積特性	記載なし	先端部の堆積が後続流の移動性を低下させる
表面形状指標	記載なし	勾配・粗度の平均と標準偏差が流動特性を反映

技術的課題と今後の展望

項目	土木研究所資料	Geomorphology論文
データ取得の安定性	・HDD→SSD換装による改善提案 ・温度上昇対策 ・リモート監視の必要性	・霧・豪雨時の検出範囲制限 ・3m以内の点群歪み ・観測継続性の確保
流速計測の課題	・1点のみの計測限界 ・3D-LiDARへの発展可能性	・0.1秒間隔での粒子追跡の限界 ・粒径との相関の課題
維持管理	・光学窓の火山灰付着対策 ・定期的な清掃(約2ヶ月毎)	・ワイヤーセンサーの設置・維持 ・ソーラーパネル+バッテリーの管理
今後の発展	・ROSを用いた制御プログラム開発 ・降雨連動起動システム	・自然渓流での更なるデータ蓄積 ・複数地点での同時観測

文献： Debris Flow + LiDAR +PIV+DNN

High‐Frequency 3D LiDAR Measurements of a Debris Flow: A Novel Method to Investigate the Dynamics of Full‐Scale Events in the Field - Aaron - 2023 - Geophysical Research Letters - Wiley Online Library

https://phreeqc.blogspot.com/2026/01/monitoring-debris-flow-dynamics.htmlの引用文献です。

2025発表となっていますが、YOLO部分がこちらで検討されています。
Deep-Learning-Based Object Detection and Tracking of Debris Flows in 3-D Through LiDAR-Camera Fusion | IEEE Journals & Magazine | IEEE Xplore

AI要約

背景
土石流の破壊性は、粗粒成分(巨礫が豊富な前面)と後続の液状化スラリーの相互作用に強く影響されるが、この相互作用の理解が不十分

過去の研究の限界:
現地での高品質データが不足
前面速度と表面速度の関係が実験室スケールでは観察されているが、現地スケールでは未確認
垂直速度分布の形状と進化に関する理解が乏しい
研究の必要性:数値モデルや物理モデルを制約するための高時空間分解能の現地データが必要

手法
観測地点
場所:スイス・ヴァレー州のIllgraben土石流観測ステーション
既存設備:ジオフォンアレイ、レーダー、超音波・レーザー機器、力板、複数のビデオカメラ観測イベント:2021年9月19日に発生した約30分間の土石流
LiDAR計測システム
使用機種:Ouster OS1-64 Gen. 1
設置位置:河道中央の堰堤上方6 mに設置
仕様:
視野角:33.5°
空間分解能:64スキャンライン
時間分解能:10 Hz(1秒間に10スキャン)
1ラインあたり2,048点
トリガー:上流のジオフォンによる起動
電源の種類:記載なし

データ処理手法
特徴検出と追跡
Matlabツールボックス"groundTruthLabeler"を使用して3Dバウンディングボックスで手動ラベリング

自動速度推定: 2つの方法を開発
ヒルシェード法:点群をヒルシェード投影し、PIV(粒子画像流速測定法)で2D速度場を導出後、3D速度に投影
LiDAR-カメラ融合法: ビデオデータにPIVを適用し、LiDARとカメラの変換を推定して3D速度場を取得
全自動速度は2秒移動平均で平滑化

機械学習による粒子検出: YOLO-v5をベースとした畳み込みニューラルネットワーク(CNN)を訓練し、ビデオ映像から礫と木質残骸を自動検出(精度0.8、再現率0.6)

流動深さと前面軌跡
イベント前の河床標高とLiDAR計測データの差分から垂直流動深さを推定
前面位置を約2秒ごとに特定し、前面および後方(2 m、6 m、10 m)の表面速度を抽出

結果
先端部と表面速度の関係
先端速度は砂防ダム付近で約0.8 m/s、上流で約2 m/sに変化
先端後方6〜10 mの表面速度は先端速度の約1.5〜2倍で、先端に近づくと減速
先端が砂防ダムに接近して減速する区間(01:53〜01:55分、センサーの6〜7 m上流)では、横断方向の速度勾配が発生し側堤が形成・再移動

流動深度と速度の時間変化
先端通過後、流動深度は最大1.5 mに達し、その後速度と礫検出数とともに徐々に減少
y = 16 mからy = 5 m(砂防ダム方向)にかけて流動深度が顕著に減少(水理学的ドローダウン効果)
イベント開始約7分後に第2のサージが到達し、流動深度と速度が15秒間増加(「速度ジャンプ」)

個別粒子の運動
29個の地物(木質残骸9個、転がる礫20個)を手動で追跡
速度ジャンプ前: 転がる礫と木質残骸はほぼ同速度で移動
速度ジャンプ後: 転がる礫の速度は木質残骸の約0.6〜0.7倍に低下
先端部の礫は30 m区間全体で追跡可能で、再循環していないことが示唆された

考察
方法論的意義
空間・時間分解能:従来研究より1桁以上高い分解能を達成
新しい視点:前面維持・伝播、サージ発達、垂直速度分布に関する新たな知見

前面維持メカニズム
表面速度>前面速度:後方の表面速度が前面速度の1.5倍で、巨礫が前面に優先的に移動
巨礫の挙動:
前面に接近すると減速し、到達後は前面の一部となるか堤防を形成
再循環は発生しない(従来の実験・数値研究と異なる)

メカニズムの解釈:
前面の巨礫サイズが流動深さと同程度
前面が「ふるい」として機能し、水と細粒粒子を逃がす
間隙水圧の排水により、巨礫が前面に到達すると速度が低下
後続物質に押されて転がりと滑りの複合運動

垂直速度分布の時間変化
観測の解釈:異なるサイズの特徴(コブル、巨礫、樹木)が異なる深さまで達し、垂直速度分布の異なる位置をサンプリング

速度ジャンプ前:ブロック滑り型速度分布
速度ジャンプ後:内部せん断を伴う速度分布(単純せん断とブロック滑りの中間)

転がる巨礫と木質デブリの速度比0.6~0.7
遷移の特徴:15~30秒の間に急激に発生

制御要因:
粗粒粒子の濃度(流動深さに近いサイズの粗粒粒子は垂直せん断を抑制)
含水量の変化
水理学的効果
堰堤上流の流動深さ減少:約15 m上流まで水理的引き下げ効果が影響

今後の考慮事項:フルード数依存性

限界と今後の研究
本研究の限界:単一イベントの観測
今後の必要性:同じ観測システムでの追加イベントデータ収集による一般性の確認

巨礫と木などの位置と大きさの違いから、深度方向の速度プロファイルを推定する工夫が新規性でしょうか。LiDARの普及が現象を正しく理解する一歩につながるのかもしれません。

文献： Debris Flow + LiDAR

Monitoring debris flow dynamics: insights from 4D-LiDAR observations in Ohya landslide, Central Japan - ScienceDirect

土石流の形状をLiDARで取得する研究です。国内では土研さんがまとめていらっしゃいましたね。

AI要約

背景・研究目的
土石流は破壊力、高速度、長距離流下を特徴とし、重大な災害をもたらす。従来の観測手法(水圧センサー、ビデオカメラ、ステージセンサー、地震計など)では、土石流の急速な時空間変化を十分に捉えられなかった。
近年、自動車用LiDAR技術の進歩により、4D-LiDAR(3次元+時間)観測が可能となった。スイスのIllgrabenなどで先行研究があるが、主に砂防堰堤間での観測に限られており、自然河道での流動特性は未解明であった。

研究目的:
自然河道を流下する土石流サージの4次元形態と流動性の解明
流動表面の形態分析による流動状態(巨礫の混入、層流・乱流状態)評価の適用可能性の検討
土石流停止段階における動態の解明

調査地
場所: 日本中部の大谷崩れ地すべり跡地、一ノ沢流域
標高: 1205〜1905 m
流域面積: 0.3 km²、全流路長1000 m
年間降水量: 約3400 mm
土石流発生頻度: 年平均5回(2016〜2023年)
地質: 中新世の砂岩と破砕された頁岩の互層

観測機器・手法
4D-LiDARシステム
土石流が到達すると自動的に観測が開始され、0.1秒間隔で連続的に3次元データを記録するシステム
LiDARセンサー: Livox Horizon (Livox Corp.製)
設置地点: 上流地点(Site U)と下流地点(Site D)の2箇所
設置時期: 2023年5月
起動システム: ワイヤーセンサー + プログラマブルリレー(Omron ZEN)
制御: オンボードコンピューター(Raspberry Pi)
電源: バッテリー + ソーラーパネル + チャージコントローラー
記録時間: 土石流到達後最低2日間の連続観測
データ保存間隔: 1時間ごと
点群密度: 0.0375 pts. cm⁻²

ビデオカメラ
機種: Sony HDR-CX470
レームレート: 60 fps
解像度: 1920×1080ピクセル
設置: Site UとSite Dの両地点
用途: LiDARデータの検証、粒径測定、流動状態の判別

UAV測量
機種: DJI Phantom 4 RTK
飛行高度: 地上50〜100 m
点群間隔: 0.02〜0.10 m
DEMグリッドサイズ: 0.1 m
ソフトウェア: Agisoft Metashape(SfM解析)
GNSS: Hemisphere A52、R320

データ解析
ソフトウェア: CloudCompare (version 2.13)、QGIS (version 3.16)
分析対象: 0.1秒間の点群データから縦断・横断図、DEM作成
分析範囲: Site D: 2.0 m×2.0 m、Site U: 1.5 m×1.5 m

表面形態分析
異なるグリッドサイズのDEMから以下を算出:
勾配: Horn (1981)の式を使用
粗度: 3×3グリッド内の最大標高差(Wilson et al., 2007)
統計値: 平均値と標準偏差

結果
観測された土石流イベント
2023年8月3日の土石流
降雨: 総雨量66 mm、最大10分間雨量11 mm(66 mm/h相当)
観測地点: Site Uのみ(下流まで到達せず)
サージ数: 8回(16:37〜16:50)、うち6回をLiDARで捕捉
流動深: 約0.5〜2 mで急激に変動
流動タイプ: 全てのサージで部分飽和流の後に完全飽和流が続いた
河床変動: 16:38に堆積、16:40と16:43に侵食を確認

2023年8月14日の土石流
降雨: 総雨量242 mm、最大10分間雨量11 mm(66 mm/h相当)
Site U: 21:00〜21:07に2回、21:30〜21:45に8回のサージ、その後2 m以上の侵食
Site D: 21:41に明確なサージと堆積を検出
UAVデータ: 上流の発生域で3 m以上の侵食を確認
Site U周辺: 上流・下流約100 mにわたり3 m以上の侵食
Site D周辺: 上流30 m〜下流100 mにわたり3 m以上の堆積

流動形態
縦断・横断形状
8月3日、Site U、16:39の最長サージ:
先端部(1〜2秒): 横断面が凸型、縦断面は土石流ローブ状、部分飽和流
後続部(3〜5秒): 横断・縦断とも滑らかなプロファイル
凸型の出現: 6回のサージ中3回で確認
凸型から滑らかへの変化: 流側への堆積と関連

8月14日、Site D:
凸型の横断面は観測されず
最大サージの先端部でスプラッシュを確認
完全飽和流または高濃度流が卓越

表面形態の定量分析
勾配の平均値
先端・中部: グリッドサイズ0.06 mまで減少傾向、以降安定
後部: グリッドサイズに関わらずほぼ一定、先端・中部よりやや小さい

ビデオカメラから得た平均粒径(0.15〜0.30 m)と標準偏差の収束値がほぼ一致(後部を除く)
粗度の平均値: 全区間でグリッドサイズに比例して増加
粗度の標準偏差: 先端・中部: 全体的に後部より小さい値
グリッドサイズ<0.10 m: 中部>後部>先端
グリッドサイズ>0.10 m: 後部で増加傾向、先端・中部は増加せず

粒径と形態指標の相関
積分値と粒径: p=0.06(有意ではないが傾向あり)
積分値/最大グリッドサイズと粒径: p=0.02(有意な相関)
グリッドサイズ0.10 mの標準偏差と粒径: 相関なし(p>0.10)

堆積過程
8月3日、Site U、16:38:
縦断方向約12 mにわたり堆積
逆勾配を形成し、10秒間でバックステッピング
約1 mの堆積高
先端部が土石流到達前の縦断プロファイルとほぼ平行

8月14日、Site D、21:41:
縦断方向約8 mにわたり堆積
逆勾配を形成し、20秒間でバックステッピング
堆積先端部が到達前のプロファイルとほぼ平行、後部で逆勾配が顕著
堆積速度: Site Uでは縦断方向に一定、Site Dでは4〜8秒後に大量堆積、その後減少

考察
流動特性
横断面形状の違い
部分飽和流の凸型形状: 河岸との摩擦影響が少ない中央部での選択的流動が原因
完全飽和流での凸型の不在: Site Dでスプラッシュのみ観測
内部力の違い: 部分飽和流は摩擦力が支配的(Oya et al., 2024)、凸型形成には摩擦力の優位性が必要

表面形態指標の解釈
勾配・粗度の標準偏差が後部で高い理由: 飽和流の乱流による表面擾乱ビデオ画像との比較: 先端部は巨礫で覆われた層流、後部は泥水主体の乱流
粒径の影響: 部分飽和流で巨礫に完全に覆われている場合、粒径も勾配の標準偏差に影響
点群密度の制約: 0.0375 pts. cm⁻²は不十分、Wang et al. (2013)は直径>63 mmの巨礫測定に1 pts. cm⁻²が必要と報告
0.1秒間隔での巨礫移動: ビデオ画像から求めた粒径との相関が低くなった原因

グリッドサイズの影響
小さいグリッドサイズ: 先端・中部・後部の違いが不明瞭
乱流の捕捉: 十分大きいグリッドサイズが必要
粒径分析: グリッドサイズが粒径より大幅に大きいと個々の粒子の影響が不明瞭
結論: 研究対象の特性に応じた適切なグリッドサイズの選択が必要

イベント特性
流下距離の違い
8月3日: 下流Site Dに到達せず、短い降雨(2.3時間)、部分飽和流が卓越
8月14日: 下流まで到達、長い降雨、完全飽和流が卓越

流動性の違い: 部分飽和流は間隙が完全に流体で満たされていないため、粒子間および流れと河道間の摩擦が増加し、流動性が低い(Major, 1997, 2000; Major and Iverson, 1999)

8月14日の高い流動性: 上流での2 m以上の侵食により大量の土砂供給、河道勾配の減少に伴う流動性低下を補完

堆積特性
逆勾配地形の形成: 堆積先端部が元の河床に対して逆勾配を形成
バックステッピング: 逆勾配が後続流の流動性を低下させ、河道での土砂の後方充填をもたらす
堆積速度の違い: 土石流の内部力の違いを反映、今後の物理メカニズム議論に活用可能

結論と今後の課題
成果
4D-LiDAR観測により以下が明らかになった:
流動高の空間的変化は土石流表面形態より河床形態に影響される
完全飽和流と部分飽和流で縦断・横断形状、表面形態、堆積特性が異なる
表面形態(勾配・粗度の標準偏差)から乱流状態や粒径の評価が可能
堆積先端部の逆勾配地形が後続流の流動性を低下させ、バックステッピングを引き起こす

限界と課題
時間解像度: 0.1秒間隔では急速な変化を見逃す可能性
点群密度: 0.0375 pts. cm⁻²では粒径分布評価の精度に限界
補完データの必要性: 地盤振動計、地震計、光ファイバーなどとの統合が望ましい
高密度・短間隔データの必要性: より詳細な粒径分布と表面形態の解明に必要

応用可能性
警報システムでの土石流検知
河川管理における河床レベル監視
災害軽減への貢献

ラズパイでも消費電力はそこそこありますし、LiDARを連続で動かして、となるとかなりの電力が必要でしょう。ソーラーパネルとバッテリーは大きなものになりそうです。林の中では難しそうです。2日毎にバッテリー交換という計画も現実的ないでしょうから、連続観測を実施したい場合は商用電源を引くしかないでしょうか。

水理公式集例題集2024

 これもプロに教えていただいたのですが、水理公式集例題集のプログラムがGitHubで公開されていました。1年間、気づきませんでした。

土木学会2024年12月新刊のご案内『水理公式集例題集（2024年版）』 | JSCE.jp for Engineers

地下水流動はFortranのままでしたが、不定流は Python になっていました。時代でしょうか。

不定流、厳しい条件だと振動していましたが、コンパイルせずに試せるようになったのはGood！お手軽に試すことができるようになりました。

摩擦項の半陰解法

前回の続きです。

1次元なのに発散することがあり、試行錯誤した結果、下記で安定しました。1次元だと思って舐めていました。

空間離散化

• 連続式: Upwind法（流れ方向に応じた風上差分）
• 運動量式: 中心差分

時間積分
• 連続式: 陽解法
• 運動量式: 半陰解法（Semi-implicit法）

効いたのは摩擦項への半陰解放の適用。摩擦項Sfの非線形性により陽解法では不安定になりやすいため、以下の線形化を実施しました。

Manning式のSfを現在の流速u_iの周りでTaylor展開：
Sf(u) ≈ Sf(u_i) + (u - u_i)·∂Sf/∂u|_{u_i}
     = Sf_const + Sf_coef·(u - u_i)

ここで：
Sf_const = n²·u_i·|u_i| / Rh^(4/3)          （定数項）
Sf_coef  = 2·n²·|u_i| / Rh^(4/3)            （線形係数）

運動量式の離散化
元の式：
∂u/∂t + u·∂u/∂x = -g·∂h/∂x + g(S₀ - Sf)
離散化：
(u_new - u_i)/dt + u_i·du_dx = -g·dh_dx + g·S₀ - g·Sf

線形化した摩擦項を代入して整理すると：
u_new·(1 + dt·g·Sf_coef) = u_i - dt·(u_i·du_dx + g·dh_dx - g·S₀ + g·Sf_const - g·sf_coef·u_i)
u_new = RHS / (1 + dt·g·Sf_coef)
RHS = u_i - dt·(u_i·du_dx + g·dh_dx - g·S₀ + g·Sf_const - g·sf_coef·u_i)

あとでプロと話していたら、1次元で不安定になりやすいのはやはり摩擦項で、半陰解法を用いて安定化させるのは実務でよくやる対応らしいです。1次元だから皆さん同じ方向に向かうのかもしれません。

1次元浅水流方程式の表現

1次元浅水流方程式（Saint-Venant 方程式） を用いて、矩形断面開水路における非定常流れを計算しました。上流端から一定時間流量を与えた後、流入を停止した際の流量減衰曲線（Recession Curve） をいくつかの条件で見てみたかったのです。

連続式と運動量式を組みあわせるだけなのですが、いくつかの表現があります。

連続式（質量保存則）
∂A/∂t + ∂Q/∂x = 0
∂h/∂t + ∂(uh)/∂x = 0
∂h/∂t + (1/B)·∂Q/∂x = 0

ここで、 h : 水深 [m]、u : 断面平均流速 [m/s]、Q : 流量 [m³/s]、A : 流水断面積 [m²] (= B·h)、B : 水路幅 [m]、t : 時間 [s]、x : 流下方向距離 [m]

 運動量式（運動量保存則）
∂Q/∂t + ∂(Q²/A)/∂x = -gA·∂H/∂x + gASf
∂Q/∂t + ∂(Q²/A)/∂x = -gA·∂h/∂x + gA(S₀ - Sf)
∂u/∂t + u·∂u/∂x = -g·∂h/∂x + g(S₀ - Sf)

ここで、g : 重力加速度 = 9.81 m/s²、H : h+z [m]、S₀ : 河床勾配（正：下流に向けて下がる）[-]、 Sf : 摩擦勾配（エネルギー線勾配）[-]

摩擦勾配 Sf は Manningの粗度式で表されます。
Sf = (n²·|u|·u) / R^(4/3)
Sf = (n²·Q·|Q|) / (A²·R^(4/3))

ここで、n : Manningの粗度係数 [s/m^(1/3)]、R : 径深（潤辺）[m]

矩形断面の径深
R = A/P = (B·h)/(B + 2h)
幅広水路（B >> h）の場合
R ≈ h

数学ではいずれも等価なのですが、時間方向に陽解法を用いると微妙に結果が異なることもあります。一番計算手順が少なくなる（誤差の積み重ねが少なくなる）表現を選択するのが良いのでしょう。

集水井の目詰り

土研さんから集水井の目詰りに関する資料が出ていました。

茶臼山地すべりの集水井における集水ボーリングの目詰まりの実態調査-論文・刊行物検索－土木研究所

2年間にわたり、目詰りの様子を写真に撮られていました。反応速度論も議論できるかなと思いつつ見ていたのですが、残念ながら一気に詰まっていました。
水質はグラフと最大、最小、平均値で示されています。どの最大値かなと迷いましたが、平均値があるので、おそらく時期はバラバラで各項目の最大値なのでしょう。現在の地すべりの専門家は水質に興味がないのかもしれません。
グラフを見る限り水質はあまり変化していなさそうでしたので、昔の文献を探れば平均的な水質は分かりそうでした。が、文献を孫引きしても水質は載っていません。おそらく、引用を間違われたのでしょう。多分、これです。残念ながらこちらにはFeが載っていませんでした。
茶臼山地すべり地上部の地下水について(II)

目詰りの有無については、沈殿量を計算すればある程度目安が付くと考えています。以下は2010年ですから、もう15年前。これにも書いていますが、地すべりの場合は詰まるとわかってもそこに打つしかないので、進展がないのも仕方がないのでしょう。
https://phreeqc.blogspot.com/2010/12/2.html
この１つは土研の方が書かれていましたが、水質の情報が揃っています。昔の地すべりの研究者は基礎を疎かにされていないようで、見習わねばなりません。

この集水井の目詰りの研究？、今後どのような着地を見せるのでしょうか。

地すべり画像と機械学習

画像を与えて機械に説明させるという手法は昨年から導入しています。

一昨年まで画像ならCNN！とか言っていたのに、ViTにその座を奪われつつあり、移り変わりが速いなあと思っていたら、昨年末には一般誌にコンクリートのクラック画像をLLMで説明させる内容が掲載されていました。もう、世間はそこまで進んでいるのだなとプチショックを受けていたところです。ちょっと気を抜くと取り残されそうで、怖いですね。

Computer Vision の専門家からも、提案がありました。
マルチモーダルAIによる地すべり画像解析と災害リスク評価 – Tohoku CVLab

Fig13 を見ると、全体としてはまだ追いつかれていないかな、と少しばかり安堵しますが、時間の問題でしょう。右側のVQA-LLMはなかなか。GTのLandslideより、土石流の方が正解のように見えます。左側は天然ダムと間違えるのも感覚的にわかります。面白い。

専門家でなくとも専門家のような答えを出せる機械、これは専門家も欲しいのではないでしょうか？近い将来における機械学習結果の使い方が楽しみになります。

Synthetic Control

Synthetic Control は、ベイズ統計を用いた因果推論です。機械学習ではなく、MCMCです。

1. モデル設定（合成コントロールの基本形）

ある「介入クラス」A について、時点 t の平均テスト点数を y_A,t とする。また、介入を受けていない対照クラスを j = 1,…,J とし、それぞれの平均点を y_j,t と書く。
合成コントロールの考え方は「A クラスは、対照クラスの“重み付き平均”で近似できる」とみなすことである。そこで、時点 t における A クラスの期待される点数（＝介入がなかったときの基準）を

μ_t = w_1 y_1,t + … + w_J y_J,t = Σ_j w_j y_j,t

と表す。w_j は 0 以上で和が 1 になるような重みとする。
実際の点数 y_A,t は、体調・偶然・測定誤差などにより μ_t からずれる。このズレを「平均 0、標準偏差 σ の正規分布に従う誤差」と仮定して

y_A,t = μ_t + 誤差_t
誤差_t ~ Normal(0, σ)

とおく。確率変数として書けば

y_A,t|w, σ ~ Normal(μ_t, σ)

である。

2. ベイズの定理と記号

推定したいパラメータを

θ = (w_1,…,w_J, σ)

とまとめて書く。ベイズの定理は以下となる。

p(θ|データ) = p(データ|θ) p(θ) / p(データ)

p(θ|データ)：事後分布（データを見た後に、θ がどの程度になりそうか）
p(データ|θ)：尤度（θ が真なら、このデータはどれくらい出やすいか）
p(θ)：事前分布（データを見る前に、θ がどのような値を取りやすいと考えるか）
p(データ)：正規化定数（全部の θ について積分した値）

実際の計算では p(データ) を明示的に計算する必要はなく、

p(θ|データ) ∝ p(データ|θ) p(θ)

という比例関係として扱うことが多い。

3. 事前分布 p(θ)

重みベクトル w = (w_1,…,w_J) について

w ~ Dirichlet(α_1,…,α_J)

とおく。誤差のばらつき σ については

σ ~ HalfNormal(τ)

とする。これらをまとめると、事前分布は

p(θ) = p(w_1,…,w_J, σ)
= Dirichlet(w | α_1,…,α_J) × HalfNormal(σ | τ)

となる。ここでは w と σ を互いに独立と仮定している。

4. 尤度 p(データ|θ)

合成コントロールの重み w は、「介入前」で学習する。介入前の時点を t = 1,…,T0 とし、A クラスの平均点：y_A,t、各対照クラスの平均点：y_j,t が観測されているとする。モデルでは、各 t について

y_A,t|θ ~ Normal(μ_t, σ)
μ_t = Σ_j w_j y_j,t

であるから、ある θ が与えられたとき、「その θ のもとで、実際に観測された系列 {y_A,1,…,y_A,T0} がどれくらい出やすいか」を表すのが尤度である。
時点 t ごとに独立と仮定すると、介入前データ全体に対する尤度は

p(データ|θ) = Π_{t=1}^{T0} p(y_A,t|θ)
= Π_{t=1}^{T0} Normal(y_A,t ; μ_t, σ)

となる。ここで μ_t は w と対照クラスの観測値 {y_j,t} から一意に決まる。
・μ_t に近い y_A,t が多ければ、この θ の尤度は大きくなる
・μ_t から遠く離れた y_A,t が多ければ、尤度は小さくなる
という直感と対応している。

5. 事後分布 p(θ|データ)

ベイズの定理より、介入前データを見た後の θ の分布（事後分布）は

p(θ | データ) ∝ p(データ | θ) p(θ)
= {Π_{t=1}^{T0} Normal(y_A,t ; μ_t, σ)}
　× Dirichlet(w | α_1,…,α_J) × HalfNormal(σ | τ)

で与えられる。これは介入前の A クラスの点数系列をどれだけよく再現するか（尤度）、もともと重み w や σ をどう想定していたか（事前）の両方をバランスさせた結果である。
この事後分布を式のまま解析的に求めることは難しいため、MCMC（たとえば NUTS というアルゴリズム）を用いて

θ^(1), θ^(2), …, θ^(S) = (w_1^(s),…,w_J^(s), σ^(s)), s=1,…,S

という多数のサンプルを生成し、p(θ|データ) を近似する。それぞれの θ^(s) が「あり得る合成コントロールの重みと誤差の大きさ」を表しており、全体として「どのような w がどの程度ありそうか」を数値的に把握できる。

6. Synthetic Control による介入効果推定

6.1 介入「前」のデータから事後分布 p(θ|データ) を推定
上で述べたように、介入前 t=1,…,T0 のデータを用いて事後分布 p(θ | データ) を求める。

6.2 介入「後」のシンセティック点数 μ_t^(s) を計算
介入後 t = T0+1,…,T の期間に注目する。この期間についても、各対照クラス j の平均点 y_j,t が観測されているとする。それぞれの事後サンプル θ^(s) に対して、時点 t ごとに

μ_post,t^(s) = Σ_j w_j^(s) y_j,t

を計算する。これは「サンプル s に対応する合成コントロールの重みを用いて、介入後の A クラスの“カウンターファクチュアル（介入がなかった場合の想定点数）”を計算したもの」である。したがって、{μ_post,t^(s)}_{s=1,…,S} は「介入がなかったとしたとき、時点 t における A クラスの点数はどのくらいになり得るか」を表す事後予測分布となる。

6.3 介入効果の事後分布を得る
実際に観測された介入後の A クラスの平均点を y_A,t（t > T0）とすると、サンプル s における時点 t の介入効果は

effect_t^(s) = y_A,t − μ_post,t^(s)

と定義できる。これは
・正なら「実際の点数の方が高い（介入がプラス効果）」
・負なら「実際の点数の方が低い（介入がマイナス効果）」
ことを意味する。

時点 t ごとに {effect_t^(s)} の平均
 → その時点の平均介入効果 2.5%点と 97.5%点
 → 95% 信用区間（事後分布の幅）を
求めれば、「いつ」「どれくらい」介入効果が表れたかをベイズ的に評価できる。
さらに、介入後の全期間について effect_t^(s) を平均すれば
・介入後期間全体の平均効果
・その信頼区間（信用区間）
も同じ要領で推定できる。

このようにして、従来の Synthetic Control（複数の対照ユニットの重み付き平均）を「事前分布」「尤度」「事後分布」の枠組みの中にきちんと位置づけ、ベイズ的な不確実性評価と一体化した形で介入効果を推定する。

7. Pythonコード（AI提案）

import pandas as pd
import numpy as np
import pymc as pm
import arviz as az

# --- データ準備 ---
# yA:  shape=(T,)       介入クラスA の時系列
# Y_donors: shape=(T,J) 対照クラス群の時系列
# T0: 介入前の最後の時点（0始まりのインデックス）

T, J = Y_donors.shape
yA_pre       = yA[: T0+1]            # 介入前データ
Y_donor_pre  = Y_donors[: T0+1, :]
yA_post      = yA[T0+1 :]            # 介入後データ
Y_donor_post = Y_donors[T0+1 :, :]

# --- 6.1 介入前データでモデルフィッティング ---
with pm.Model() as sc_model:
    # 事前分布
    w     = pm.Dirichlet("w",     a=np.ones(J))   # w_j ≥ 0, Σw_j=1
    sigma = pm.HalfNormal("sigma", sigma=1.0)     # σ>0

    # 合成コントロールの期待値 μ_t
    mu_lin = pm.math.dot(Y_donor_pre, w)
    mu_pre = pm.Deterministic("mu_pre", mu_lin)

    # 尤度
    pm.Normal("y_obs", mu=mu_pre, sigma=sigma,
              observed=yA_pre)

    # サンプリング
    trace = pm.sample(
        draws=1000,          # 事後サンプル 1000
        tune=1000,           # チューニング 1000
        chains=4,            # チェーン数
        cores=4,             # 並列数
        target_accept=0.9,   # NUTS の受諾率
        max_treedepth=12,
        random_seed=0,
        progressbar=True
    )

# --- 6.2 介入後のカウンターファクチュアルを計算 ---
# ポスターサンプルから w を (S,J) にフラット化
w_samples = trace.posterior["w"] \
    .stack(sample=("chain", "draw")) \
    .values  # shape = (S, J)

# μ_post[s,t] = Σ_j w_samples[s,j] * Y_donor_post[t,j]
# 結果 shape = (S, T_post)
mu_post = (Y_donor_post @ w_samples.T).T

# --- 6.3 介入効果を計算 ---
# effect[s,t] = 実測 yA_post[t] − μ_post[s,t]
effect = yA_post[None, :] - mu_post  # shape = (S, T_post)

# 時点別平均効果と95%信用区間
effect_mean_t = effect.mean(axis=0)                # (T_post,)
effect_hdi_t  = az.hdi(effect, hdi_prob=0.95)      # (T_post, 2)

# 介入後期間全体の平均効果とその信用区間
avg_effect     = effect.mean()                     # scalar
avg_effect_hdi = az.hdi(effect.flatten(), hdi_prob=0.95)  # (2,)

# --- 出力例 ---
print("時点別平均効果    :", effect_mean_t)
print("時点別95%信用区間:", effect_hdi_t)
print("全期間平均効果    :", avg_effect)
print("全期間95%信用区間:", avg_effect_hdi)

CausalBERT

因果推論にテキストを扱う時代になっていたのはプチ衝撃でした。

BERTの拡張なので、Transformer, PyTorchの一般的な環境ですぐに使えます。

発表された文献と読んだ図書では利用方法が異なりましたが、どちらも使えそうです。以下は後者の復習です。AIが提案するコードは私の環境ではやや古く、動作未確認です。

1. 目的と構造

目的：テキストと構造化データから、ある処置の因果効果（ATE/ITE）を推定する。

同時に2つのモデルを学習する二重ロバスト型：

処置モデル（gモデル）：g(T | X_text, X_cov)

アウトカムモデル（Qモデル）：Q(Y | T, X_text, X_cov)

ここで

X_text: テキスト特徴

X_cov: 交絡因子

T: 処置変数（0/1 など）

Y: アウトカム（結果）

2. 損失の重み g_weight と Q_weight

g_weight

処置モデルの損失重み。処置の割り当てメカニズムをどれだけ重視するかを調整。

Q_weight

アウトカムモデルの損失重み。因果効果推定の中心となるアウトカム予測をどれだけ重視するかを調整。

概念的な損失式：

Loss_total = g_weight * Loss_g + Q_weight * Loss_Q + ...

3. ATE / ITE の計算

推論結果 preds を

preds[:, 0] ≈ Ê[Y(0) | X]

preds[:, 1] ≈ Ê[Y(1) | X]

という反事実アウトカムの推定値と解釈できる場合：

個別因果効果（ITE）：

ITE_i = preds[i, 1] - preds[i, 0]

平均処置効果（ATE）：

ATE = mean(preds[:, 1] - preds[:, 0])

無交絡などの識別仮定が成立しているとき、これを

E[Y | do(T=1)] − E[Y | do(T=0)]

として解釈できる。

4. 典型的なデータ構造

text_main: 主テキスト

text_aux: 補助テキスト（例：カテゴリ名など）

treatment: 処置変数（例：A/B条件）

covariate: 交絡因子（例：属性フラグ）

outcome: アウトカム（例：クリック有無）

CausalBERT に渡すテキストは 1 本なので、複数テキストを結合して使う：

5. 学習〜推論コード（AIの提案）

import pandas as pd
from causalnlp.core.causalbert import CausalBertModel

df=pd.read_csv("data.csv")

#複数の文字列項目を 1 つのテキスト入力として利用。
df["text"] = df["text_main"].fillna("")
             + " [SEP] "
             + df["text_aux"].fillna("")

cb = CausalBertModel(
    batch_size=32,
    g_weight=0.01,
    Q_weight=1.0,
    mlm_weight=0.1,
	max_length=128
)

#学習
cb.train(
    texts=df["text"], 
    confounds=df["covariate"],
    treatments=df["treatment"],
    outcomes=df["outcome"],
    epochs=6,
)

#ATE の算出
print(cb.estimate_ate(texts=df["text"]), confounds=df["covariate"])

文献：Benchmark of FLAC3D. Comparison with 3DEC

Benchmark of FLAC3D. Comparison with 3DEC and other numerical simulation tools – SKB.com

AI要約

背景
スウェーデンの核燃料・廃棄物管理会社（SKB）は、処分場の設計と長期的な安全評価において、岩盤内で発生する機械的、熱機械的、動的プロセスを数値シミュレーションすることが不可欠であると考えている。長年にわたり、3DECがこれらの岩盤力学シミュレーションの主要なツールとして使用されてきた。しかし、モデリング作業の発展とスウェーデン放射線安全庁（SSM）からの要求により、SKBは代替の数値シミュレーションツールの使用を検討している。
FLAC3Dは、その計算エンジンの近代性と速度から、3DECの代替として検討されている。ただし、FLAC3Dにおける不連続面の取り扱いロジックは、少なくともバージョン7においては、3DECほど堅牢ではない。この研究では、FLAC3Dバージョン7と、より堅牢な新しい不連続面ロジックを搭載したベータ版のバージョン9の両方を比較する。

手法
複数のモデリングケースを検討し、FLAC3Dの結果を3DECおよび他のコードからの対応する結果と比較した。検討ケースは以下のとおり。:
単一亀裂モデル（準静的）
複数亀裂モデル（熱機械的、準静的）
開口部と亀裂を持つ熱機械モデル（準静的）
フォルスマルク地震モデル（動的）
大規模地震ベンチマークケース（動的）
すべてのモデリングケースで、さまざまなレベルの複雑さを持つジョイント面が含まれていた。いくつかのケースでは、1つのジョイントのみが含まれ、他のケースでは、複数の交差するジョイントが掘削とともに含まれていた。
熱機械的シミュレーションでは、3DECからFLAC3D_7に温度と温度増加をインポートした。FLAC3D_9では、分析的熱ロジックが実装されたため不要であった。

結果
ジョイントの少ないモデリングケースでは、FLAC3Dの両方のバージョンが3DECとほぼ同一の結果を生成した。
複数の交差するジョイント面がある場合、FLAC3Dソリューションでインタフェースロジックを使用すると、交差点の周りに局所的な変位異常が発生する可能性あり。
FLAC3Dバージョン9でゾーンジョイントロジックを使用すると、FLAC3Dバージョン7で得られるよりも信頼性が高く、高速なソリューションが得られた。
大規模地震ベンチマークケースでは、FLAC3Dバージョン9は、動的応答のある程度の過小評価を示した。

考察
FLAC3Dバージョン9は、3DECや他のコードによって生成された対応する結果と数パーセント異なる結果を生成した。コードのアプローチの違いを考えると、これは満足できると考えられる。
FLAC3Dの計算速度は3DECよりも大幅に向上しており、少数の交差するジョイント面を含むアプリケーションでは、FLAC3Dバージョン9が3DECの魅力的な代替手段になる可能性がある。
両バージョンのFLAC3Dは、多数のジョイントを持つが交差がない問題に対して、3DECの代替手段になる可能性がある。
全体として、この研究は、FLAC3Dが特定の問題領域、特にFLAC3Dバージョン9の新しいゾーンジョイントロジックを使用した場合に、3DECの実行可能な代替手段であることを示唆している。

3DEC の代わりにFLAC3Dを利用できるか？という内容の文献です。 Swedish Radiation Safety Authority（SSM）からどのような要請があったのか気になります。
インターフェース要素は変形の取り扱いが難しそうですが、Ver9で良い結果が得られている点は覚えておきましょう。

文献：Estimating groundwater and stream water ages with CFC

Estimating groundwater and stream water ages with chlorofluorocarbons – SKB

AI要約

背景
目的: 塩素化フッ化炭素 (CFC) を用いて、スウェーデン北部のクリックラン研究流域における地下水と河川水の年代を推定する。
CFCの利用: CFCは、0-50年程度の「若い」地下水のトレーサーとして有用である。過去の研究(Kolbe et al., 2020)から、CFCは地下水の年代の深度依存性も検出できる可能性が示唆された。
先行研究: 2017年の研究で、浅い氷堆石帯水層において、地下水年代が深度とともに増加する傾向が確認された。表層地下水（地表面から2-3m）でも、既に30年前のものであった。
研究の動機: 先行研究の結果を検証するため、涵養域と湧出域の両方で地下水を採取する。また、同一地点での繰り返しサンプリングによるCFC年代の整合性を検証し、河川水の年代推定への応用可能性を探る。

手法
調査対象: クリックラン流域とデゲロ・ストルミル。（ページ2）
サンプリング:
2021年: 地下水44試料 (28井戸)、河川水2試料 (C2, C7出口)
2022年: 地下水8試料 (深井戸、河畔域井戸、泥炭地)、河川水18試料 (高流量時と低流量時)
2023年: 河川水29試料 (異なる流量条件)
CFC分析: フランス・レンヌの分析プラットフォーム "CONDATE Eau" にて、パージ・アンド・トラップ気体クロマトグラフィー法で分析。
年代推定: ピストンフローモデルを用いて、CFC濃度から地下水年代を推定。
モデルの仮定:
（地下水の）水が特定の流れに沿って移動する場合、すべての水は同じ移動時間を持つ。
水サンプル内の水は、採取された流れにそって、同じ年齢を持つ。
CFCが水力学的分散や混合の影響を受けず、一定の流れ場で地下を保守的に輸送される。
均質な厚さと均一な地下水涵養を備えた理想化された不圧帯水層内で輸送が発生する。

結果
地下水年代: 地下水年代は深度とともに増加する傾向が見られた。表層地下水でも、既に数十年経過している。
再サンプリング: 再サンプリングした地点では、地下水年代が数年古くなる傾向が見られた。
河川水年代: ピストンフローモデルによる河川水年代は、35-55年の範囲にあった。高流量時や小流域では、河川水年代が若くなる傾向が見られた。
2021年のサンプリング結果: 詳細なCFC濃度と年代推定結果がTable 3-1に示されている。 サンプリング深度と不飽和帯の厚さ、CFC年代との関係はTable 3-2に示されている。

今後の提案される作業
モデリング:
ピストンフローモデルだけでなく、より複雑な（2-3パラメータの集中モデルや2D/3D物理ベースの数値モデル）モデルをテストする。これにより、平均年齢だけでなく、水年齢の分布を推定できるようになる。
安定同位体データとの組み合わせ:
酸素18や重水素などの安定同位体データ（数日から数か月オーダーの情報を提供）とCFCデータを組み合わせることで、河川水の水年齢分布に関する新たな知見が得られる可能性がある。
降雨時のCFCサンプリング:
融雪期や渇水期だけでなく、降雨時の河川水のCFCサンプリングを行うことで、水の貯留・放出メカニズムの理解を深める。
これらの提案は、CFCデータを用いた水Age推定の理解を深め，より高度な水文モデルの構築に役立つ。 

簡単な1次元輸送でもモデル化できたという結果が良いと思います。湧出箇所だと難しいかもしれませんが、まずは概念モデルを簡単なモデルで検証できたという点で現象を理解していることになるのでしょう。

文献：Groundwater flow measurements in permanently

Groundwater flow measurements in permanently installed boreholes. Test campaign 18, 2024 – SKB.com

AI要約

背景
フォルスマルク処分候補地では、天然状態での地下水流動の把握が、サイトの水理・水化学特性や人工バリア機能評価に重要。
2005年以降、永久設置ボーリング孔で年1回の長期監視を継続。
2013–2017年の解析で、蒸発や測定時間の影響が大きいこと、旧サンプリング装置の劣化が問題と判明。
2019年以降、オンライン蛍光計測方式を導入し、2020年からは完全移行。2024年は第18回キャンペーンで25区間を測定。

手法（ダイリューション法の詳細）
原理：パッカーで隔離した区間に蛍光トレーサ溶液（Amino‑G Acid）を注入し、区間内を循環させて均一濃度（約0.5–1 ppm）にした後、自然流入地下水による濃度低下を連続計測する。
装置：区間と地表をつなぐ循環・注入用チューブと圧力管を用い、ダウンホールポンプで閉回路循環を行う。GGUN‑FL30蛍光計測器で数秒～数分間隔のオンライン測定を実施し、濃度は事前の二点校正（100・1000 ppb）に基づきmV→ppbに変換する。
試験操作：
まず15–20分循環し、オンラインで背景蛍光を確認・固定。
セクション1容積を循環するのに要する時間と同じ時間だけ、1 000 mg/l のトレーサ原液を注入（注入/循環流量比1/1000）して均一濃度を形成。
以後、5分間隔で長期間（1–3週）濃度を記録。試験終了後は少なくとも区間体積の1倍以上を揚水し残留トレーサを除去。
データ処理：
ガス気泡によるスパイクを移動平均±5 ppbから負側に外れる点として除去。
希釈相の ln c – 時間プロットから、注入直後の非定常部分を避けて直線部分を視認で選定し、回帰直線の傾きから式 ln(c/c0) = − Qbh/V · t を用い、区間容積 V を掛けて孔内流量 Qbh を算出。Qbh をボーリング径・区間長・補正係数 α（割れ目岩で一律2と仮定）で補正し、ダーシー流速 v を求め、既存の透水係数 K と組み合わせて水理勾配 I=v/K を推定する。

結果
2024年は25区間を測定し、流量0.01–31 ml/min、ダーシー流速 9.7×10⁻¹¹～1.4×10⁻⁷ m/s、水理勾配0.0001～0.7 m/m を得た。多くの区間で2020年以降の結果と整合する一方、KFM02A:3 では異常に高い流量が得られた。

考察
オンライン方式により、従来法で問題だった採水に伴う追加希釈・蒸発誤差が除去され、長期的な直線部の抽出が可能になった。ただし α の仮定や K 値の不確かさ、単一割れ目支配区間ではダーシーモデル自体の妥当性に限界があり、一部区間で異常に大きな水理勾配が推定されることが示されている。

岩盤の透水性を原位置で求める知らなかった方法です。水理勾配を小さくできるので、丁寧な測定ができそうです。

文献：montmorillonite dissolution

Reactive transport modelling of montmorillonite dissolution. Report for the safety evaluation SE-SFL – SKB.com

AI要約

背景
SFLのBHAでは、廃棄物を収めたコンクリート構造物をベントナイト（MX-80）で取り囲む設計であり、セメント系高pH間隙水との相互作用により主成分であるモンモリロナイトが溶解し得るため、長期安定性を評価する必要がある。モンモリロナイト溶解速度はpHと温度に依存し、低温ではpH依存性が弱まる。本検討は温度15 °Cの条件を対象とする。

手法
1次元反応輸送モデルをiCP（Comsol Multiphysics＋Phreeqc）で構築し、PHASTも検証用途で使用。反応は、コンクリート側は主に平衡、ベントナイト側のモンモリロナイトは速度論で扱う。溶質移行は飽和条件下の拡散支配（Fick型）。

二つの分析体系を構築:
簡略系（ベントナイトのみ）: 左境界にコンクリート間隙水を固定濃度（例：portlandite平衡）として与え、ベントナイト厚は2.0 m。右境界は閉鎖または地層水固定濃度の感度を実施。
完全系（コンクリート–ベントナイト相互作用）: セメント系（廃棄物領域とコンクリート構造体）とベントナイトを明示的に連成（ベントナイト厚2.3 m）。
化学モデル:
熱力学データベースはThermoChimie v9b。
モンモリロナイト溶解はSato-Oda式（pH・温度・平衡接近度依存）。速度のpH依存性は15 °Cでは弱く、pH 8→>13でも増加は小さい（約2倍）。
コンクリート側の反応は平衡を仮定（OPCの劣化シーケンス：アルカリ溶脱→portlandite溶解→C-S-Hの脱カルシウム化等）。
境界水（簡略系）の例:
portlandite平衡のコンクリート間隙水は15 °CでpH=12.835。ベントナイト孔隙水・地層水の代表組成も提示（Table 6-1）。
感度解析:
拡散係数、モンモリロナイト反応表面積（SSA）、二次鉱物集合、境界条件（portlandite vs C-S-H平衡水）、ベントナイト組成・交換反応、コンクリート側の拡散係数・組成（アルカリ有無、廃棄物領域の寄与）など、合計21ケース（簡略系と完全系）を実施。
反応表面積の考え方:
高密度のベントナイトでは反応面がエッジに局在し、マスキング効果により実効反応面積が従来想定より小さい可能性。SSA=800 m2/gを基準としつつ、BET近傾の30 m2/g、およびTeradaらの式に基づく0.03 m2/gの感度を実施。

結果
簡略系（左境界＝portlandite平衡高pH水）:
参照ケースでは、10万年後の溶解フロントは約1.5 m、ベントナイト全体のモンモリロナイト残存は約24 %（すなわち約76 %が溶解）。進行は時間の平方根にほぼ比例し、拡散支配。
拡散係数を1桁低下（例：1.2×10^-11 m^2/s）すると溶解深さは約0.5 m、溶解質量は約28.5 %に減少（参照比で約3倍遅い）。
右境界を地層水固定濃度にすると溶解フロントは約1.4 m（10 cm程度の抑制）。
反応表面積の感度：30–800 m2/gの範囲では総溶解量への影響は小さいが、0.03 m2/gでは溶解が顕著に低下。
二次鉱物の許容により、溶液がモンモリロナイトと遠い平衡状態に維持され、溶解が進みやすくなる。二次鉱物の総沈殿モル数が少ないほど溶解深さは小さくなる。
portlandite平衡水より穏やかな境界水（C-S-H1.2平衡）では溶解深さが約半減（参照比）。
完全系（コンクリート–ベントナイト連成）:
参照ケースでは、10万年後にベントナイト全体で約30 %が溶解（約70 %残存）。界面から約0.25 mは枯渇し、全厚で部分的影響。
界面pHは時間とともに低下し、約10.5–10.8に収束傾向（30–40 kyr）。拡散・供給が時間依存のため、フロント進行は単純な√t直線からわずかに外れる。
コンクリート側拡散係数の感度：Deを参照から3.5倍小さくすると溶解は約22 %、1桁大きくすると約40 %溶解。
反応表面積0.03 m2/gでは、10万年後も95 %以上が残存し、溶解は大幅に低減。
二次鉱物集合の変更（ゼオライト種の入替等）で最終溶解は参照比で約5 %程度低減にとどまるケースあり。
廃棄物領域を除外（閉鎖境界に置換）すると約25 %溶解に減少。初期アルカリの除去でも同程度の低減。ただし非保守的。
早期の孔隙閉塞傾向（簡略系で約300年、完全系で約2,000年）が計算上確認される（モデルは孔隙率未更新）。
解析モデルとの比較
簡略系に対するシュリンキングコアモデル（SCM）は、溶解深さが時間の平方根に比例し、pH=12.835境界で反応輸送結果と良好に整合（初期数千年を除く）。
コンクリート–ベントナイト連成の解析解（Neretnieks 2014）では、界面pHが約10.5となり、参照の反応輸送と比較してモンモリロナイト溶解深さを過小評価する傾向がある。溶解深さは、モンモリロナイト1 molの溶解に必要なOH^-モル数（f2）に反比例し、具体的にf2が小さいほど溶解が進む。

考察
簡略系は、コンクリートの拡散抵抗や有限のアルカリ供給量を考慮しないため過度に保守的。一方、完全系では界面pHが低下し、長期の溶解は緩和される。
溶解支配因子:
拡散係数（ベントナイト・コンクリート）：水飽和下では拡散支配で最重要。
反応表面積：30 m2/g以上では感度は小さいが、0.03 m2/gでは溶解が大幅に低下。
界面の扱い（固定境界 vs 連成）：結果に決定的影響。
二次鉱物集合：平衡接近度を通じて溶解速度に影響。
セメント側アルカリ・廃棄物領域の寄与：無視は非保守的。
総括：最も妥当な設定では、10万年でモンモリロナイトの20–70 %が溶解し得るが、最近の低い実効反応表面積の知見を導入すると溶解は顕著に低下し得る。

結論は非常にシンプルなのですが、たどり着くには多くの実験が必要だったのでしょう。比表面積や拡散係数をどのようにして設定するかは、まだまだ難しいようです。

因果探索

 因果推論にはDAGが必要です。それをデータから探索するのが因果探索。

 因果探索には「制約ベース」「スコアベース」「関数ベース」「連続最適化・勾配ベース」「時系列拡張」「ML・深層生成モデル」などがあります。代表的な仮定は以下の3つ。

• 因果的マルコフ条件（Causal Markov）; DAG G 内の任意のノードX、その親ノードの集合(Pa(X))が与えられたとき、自身の子孫でない変数すべてと条件付き独立であるという性質。「直接の原因が分かれば、それ以上さかのぼる必要はない」
• 信念性（Faithfulness）; 確率分布 P で観測される条件付き独立性は、DAG G のD分離から予測されるもの以外には存在しない。「この変数とこの変数は無関係に見える、という現象は、すべて因果グラフの構造によって説明できる」
• 因果十分性（Causal Sufficiency／潜在共通原因なし）; 観測された変数に、交絡因子となるような「隠れた変数（見えない要因）」は一つも存在しないという性質。「大事な変数は、全部見えている」

代表的な手法です。 

手法	タイプ	アルゴリズム	仮定、特徴	注意点
PC	制約ベース	完全無向グラフ→条件付き独立（CI）検定でエッジ削除→v-structureなどのルールで向き付け	マルコフ性・信念性・因果十分性、連続/離散可	高次元・サンプル数不足に弱い、潜在共通原因があると誤る（因果十分性が必要）
FCI/RFCI	制約ベース（拡張）	PCを拡張し、潜在変数を許す CI 検定＋向き付け規則で PAG; Partial Ancestral Graph（部分祖先グラフ）を構成	マルコフ性・信念性、潜在共通原因・選択バイアス許容	出力グラフが複雑、解釈困難になりやすい
GES/GIES	スコア＋貪欲探索	BIC などのスコアを最大にするよう、エッジ追加→削除を繰り返し DAG 探索	マルコフ性・信念性・因果十分性・スコア同値性・分解可能性、連続/離散可、GIES は介入データを利用可能な拡張版	局所最適解に陥る可能性、スコア依存性
LiNGAM	ICA＋回帰	非ガウス性・独立誤差を仮定しICA; Independent Component Analysis（独立成分分析）的に因果順序を推定→順に回帰で辺推定。	線形・非ガウス誤差、誤差独立、潜在変数なし版が基本形、主に連続変数向け	仮定が強く，ノイズがガウス寄りだと破綻
NOTEARS	連続最適化	h(W) =tr(exp(W⊙W))-d =0などの滑らかなacyclicity 制約付き損失を勾配法で最小化し、隣接行列を直接学習	線形SEM（ガウス誤差）想定、高次元対応	非線形拡張版は計算量大、局所解の可能性
GOLEM	連続最適化	尤度（例：線形ガウスSEM）＋スパース正則化＋DAGペナルテティ（対数尤度＋λ・DAG違反項）を最適化しDAG近似	線形SEM・大規模データ対応	ペナルティ法ゆえ制約は厳密でなく，λ調整が必要
DECI	深層生成モデル＋ベイズ推論	①VAE系の生成モデルで非線形SEM＋DAGを生成的に定式化 ②変分推論（またはMCMC）でposteriorを近似、基本はDAGの事後分布とSEMを学習し，介入データがあれば推定が強化される	非線形・複雑分布、観測＋介入データ対応、事前分布や制約を組み込みやすい	計算コスト大。モデル設計・ハイパーパラメータに敏感

専門知識の組み込み例

•  PC／FCI系：初期完全グラフのMaskで禁止エッジを除外。tier（層）情報や時間順を向き付け規則に追加。
•  GES／NOTEARS／GOLEM：Mask行列で必須エッジ・禁止エッジを固定。ペナルティ項の重みを変えてprior的に導入。
• LiNGAM：既知の部分順序（時間順など）で可変順序を固定。
• DECI：DAGやエッジの事前分布（スパース先行分布など）を自然に定義。

Pytorch などの慣れたツールを使えば、手軽に結果を得られるできるでしょう。が、妥当性の判断が難しい。数をこなして慣れるしかないのかな。