文献：Contact modelling approach to simulate water flow in a rock fracture

Contact modelling approach to simulate water flow in a rock fracture under normal loading – Modelling report of Task 10.2.2. Task 10 of SKB Task Force GWFTS – Validation approaches for groundwater flow and transport modelling with discrete features – SKB.com

AI要約

背景
未開口の自然岩盤亀裂に対し、法線荷重を変化させたときの流量を「ブラインド予測」する課題。
亀裂を有する岩試料で水理‐力学試験を実施後、亀裂を開口し両面を高精度スキャン、その後再装着して再試験。
目的：単一亀裂での荷重‐変形‐流れの関係と流路のチャンネリング、ならびに実務的な検証手法の検討。

手法
使用コード：COMSOL Multiphysics。
流れ：低レイノルズ数の定常流とし、Navier–Stokes式から慣性項を除いたStokes流として 3D 解析。
岩盤変形：線形弾性・等方・連続体とし、亀裂面の接触を「拡張ラグランジュ法」による接触条件で表現。

モデル群：
Analytical Model：平行平板間流れ（立方則）でスキャン精度が流量に与える影響を評価。
Parallel Plate Model：3D平行平板の数値解で解析解との整合とメッシュ条件を確認。
Small Surface Model：小さな実測亀裂面（Task 10.2.1）を用い、上面のx,y,z変位とz軸回転の影響をDSDによる感度解析。
Nearing Contact Model：上下岩体を初期ギャップから押し付ける接触モデル（1 MPa荷重）。
Departing Contact Model：初期的に重なりを与えた状態から上盤を押し上げる接触モデル（0,1,2,4,6,8 MPa）。
接触オフセット：数値的な「めり込み」を避けるため、上下面の最小間隔を強制（0.05–0.1 mm）。この仮定が開口と流量を増大。

結果
Analytical / Parallel Plate：開口 = スキャン精度0.035 mmに対応する流量は約0.1 ml/sで、実測流量と同程度。実測から逆算した等価平行平板開口はいずれもスキャン精度以下。
Small Surface Model：x,y,z変位はいずれも統計的に有意で、とくにz方向変位が流量に最も強く影響。z軸回転はほぼ無影響。
Nearing / Departing Contact：接触解析から得た変形後空隙でStokes流解析を実施。

いずれの荷重・接触アプローチでも、予測流量は実測より数桁大きい。接触オフセットやスキャン精度、初期位置誤差、補間誤差、弾性仮定を考慮しても、この差を説明できない。
von Mises応力は局所的に一軸圧縮強度を超え、現実には局所破壊や塑性変形が起こる可能性が高いことを示唆。

考察
スキャン精度0.035 mmでは、実測レベルの微小流量を信頼性高く再現するには不足。逆算開口が常に精度以下であるため、幾何データに根本的な限界。
上下亀裂面の水平方向移動（x,y）とより精密なインターロッキング、ならびに岩の塑性変形・破壊を無視しているため、実際より開口が大きく評価され、流量を過大評価している可能性が高い。
接触オフセットは数値上必要だが、同時にモデルに系統的バイアス（過大開口）を導入する。
既知の不確実性を考慮しても、モデルと実験の乖離は説明できず、「未知の支配的要因」が欠落していると結論。
接触モデル＋詳細幾何に基づく手法は、本課題のような極めて小さい流量の予測には実務的でなく、より小さい試料や高精度スキャン、あるいは別の概念モデルが必要と提言。
併せて、複数現象・複数検証を通じたモデルの「妥当性ランクR」を定義する検証評価指標を提案し、今後のモデル群の統合的評価に用いることを示唆。

0.035mmで精度が不足とは、驚きです。いえ、このような制度が出せるのも驚きなのですが。地層処分ですので細かい。
使われた 3 つの接触法（ペナルティ、拡張ラグランジュ、Nitsche）、どれも知りませんでした。拡張ラグランジュが最もオーバーラップ（食い込み）が小さいものの計算コストが高いそうです。オーバーラップを完全にはゼロにするには非常に細かいメッシュが必要で、現実的なメッシュでは 0.01–0.1 mm 程度の食い込みが残ると報告されていますが、それでも完全には抑えきれなかったとのこと。扱うスケールが全く異なります。

文献：Site scale modelling of groundwater evolution

Site scale modelling of groundwater evolution at SFR. A modelling application with iDP – SKB.com

AI要約

背景
スウェーデン核燃料・廃棄物管理会社（SKB）は、Forsmarkにある低・中レベル放射性廃棄物の最終処分場SFRの拡張を計画している。地質学的な処分の候補地の信頼性ある評価には、地下水の流動（hydrogeology）と地下水化学（hydrochemistry）を統合したモデル構築が必要である。しかし、数値的・概念的な制約や計算負荷のため、多くの従来研究ではこれらの重要なプロセスを個別に扱ってきた。本研究では、SFR周辺の地圏の進化を記述する3次元反応輸送モデルの構築と検証を目的とし、地下水流動と化学反応を統合したモデル化に取り組んだ 。

手法
本研究では、DarcyToolsとPFLOTRANという2つの計算コードを連結する「iDP（interface DarcyTools-PFLOTRAN）」を活用した。DarcyToolsは複雑な地層中の流れ・物質輸送を、PFLOTRANは高性能計算環境での流れ・反応輸送をそれぞれ扱う。DarcyToolsで得られた透水性・圧力場等をPFLOTRAN形式に変換し、反応輸送シミュレーションを行う。地圏を等価連続多孔質媒体（ECPM）として上流から下流までを模擬し、rock matrix領域は「multi-continuum approach（Lichtner法）」により反応性を持たせている。
モデルの検証のため、1D、2D、多孔質岩石行列を想定したモデルや解析解を利用。化学反応は前研究FASTREACTとの比較のためベースケース（方解石＋ヘマタイト平衡）とバリアントケース（方解石＋非晶質硫化鉄(FeS(am))平衡）の2通りを設定。3500年間（2500〜6000AD）の沿岸線移動および気候変動・地下水組成の変化も考慮した。

結果
3D反応輸送計算により、岩石マトリクス領域を含むモデルと含まないモデル（FASTREACTと同水準の設定）の比較が可能となった。
塩化物（Cl）の分布は、希釈・海岸線移動による水侵入の影響をよく再現し、岩石行列の拡散効果による急激な化学変化の緩和が見られた。
pH値は両ケースとも7.2〜8で安定。pe値、還元状態も3500年の計算期間で安定して保持された。
Fe2+等については、岩石マトリクスでの biotite 溶解による供給があり、反応的な鉱物（ヘマタイトあるいはFeS(am)）の存在条件によって、濃度推移が異なった。
主要カチオンの動態、水理パターン、時空間分布についても詳細な可視化を3Dで実現した。

Fe^2+の濃度進展には岩石マトリクスの反応性の影響が顕著であった。バリアントケース（FeS(am)平衡）ではFe^2+の濃度はベースケース（ヘマタイト平衡）に比べて低下し、FeS(am)の沈殿が効果的なFe^2+の吸収源＝シンクとなっていることが示された。 

Cl^-のような保存的成分（反応しない成分）は水収支および拡散・移流の指標として有用であり、モデルでの到達遅延や空間分布が確認された

考察
岩石マトリクス（そのパラメータ）の設定によって、地下水組成の急激な変化（例：希釈効果、流動路の変化）が拡散・反応により緩和されることが示された。
従来の1DモデルやFASTREACTによる結果と、今回の3Dモデル（岩石マトリクスを加味したもの）の比較で、岩石マトリクスを考慮することの重要性が強調された。
モデルパラメータ（特に岩石マトリクスの空隙率、表面積等）が地質環境ごとに大きく異なるため、今後は現場データによるパラメータ同定が重要になる。
開発した3Dモデル（iDP）は、今後より複雑・広範なスケールでの安全評価や、放射性核種の挙動・バリア相互作用のモデル化が期待される。

亀裂（fracture）の移流分散と、岩石マトリクス拡散を解いているようです。岩石中の拡散パラメータを設定するのは難しいでしょうね。
反応計算の地球化学的パラメータは前研究FASTREACT（Román-Ross et al. 2014）のモデルに準拠とのこと。どのようなものでしょうか。

1次元移流分散（有限差分法、Python）

GitHubで水理公式集例題集2024 のコードを眺めていた時に、地下水の不飽和浸透や移流分散計算がUPされているのを見かけました。

EXCEL版 と Fortran版 があり、懐かしいなあと思いつつ「例題 1.24 水平一次元帯水層を対象とした溶質輸送過程の計算（数値解）」の Fortran コードを Python に移植してみました。

で、結果が全く同じだったことには少し驚かされました。

以下、コードに沿ったFTCS（Forward-Time Centered-Space）による離散化とスクリプトです。手元に例題集が無いため、表現は異なるかもしれません。スクリプトではスライスを利用することで計算が1行となり、見やすくなりました。実際の地下では流速や流量が一定ではないので、ここまで簡素にはできませんが、まずはここからなのでしょう。

1. 支配方程式

∂c/∂t + v·∂c/∂x = D·∂²c/∂x² ……(1)

ここで、

c : 濃度 (mg/L)

t : 時間 (day)

x : 距離 (m)

v : 移流速度 (m/day)

D : 分散係数 (m²/day)

2. 差分近似

2.1 時間微分の前進差分

時刻 t における時間微分を前進差分で近似すると、

∂c/∂t ≈ (c_i^(t+1) - c_i^t) / Δt ……(2)

ここで、

c_i^t : 格子点 x_i、時刻 t における濃度

Δt : 時間ステップ

2.2 空間2次微分の中心差分（拡散項）

格子点 x_i における空間2次微分を中心差分で近似すると、

∂²c/∂x² ≈ (c_(i+1)^t - 2·c_i^t + c_(i-1)^t) / Δx² ……(3)

ここで、

Δx : 空間格子間隔

2.3 空間1次微分の後退差分（移流項）

v > 0（下流方向の流れ）を仮定し、上流差分（upwind差分）を適用すると、

∂c/∂x ≈ (c_i^t - c_(i-1)^t) / Δx ……(4)

3. 差分方程式の導出

式(1)に式(2)、(3)、(4)を代入すると、

(c_i^(t+1) - c_i^t) / Δt + v·(c_i^t - c_(i-1)^t) / Δx

= D·(c_(i+1)^t - 2·c_i^t + c_(i-1)^t) / Δx² ……(5)

式(5)の両辺に Δt を掛け、c_i^(t+1) について解くと、

c_i^(t+1) = c_i^t

+ (D·Δt/Δx²)·(c_(i+1)^t - 2·c_i^t + c_(i-1)^t)

- (v·Δt/Δx)·(c_i^t - c_(i-1)^t) ……(6)

ここで、

r_diff = D·Δt / Δx² （拡散数）

r_adv = v·Δt / Δx （クーラン数）

を用いると、式(6)は以下となる。

c_i^(t+1) = c_i^t

+ r_diff·(c_(i+1)^t - 2·c_i^t + c_(i-1)^t)

- r_adv·(c_i^t - c_(i-1)^t) ……(7)

import numpy as np
import csv

# +++++++++++++++++++++++++++++++++++++++++++++++
# 1次元移流分散方程式の数値解析
# 有限差分法：FTCS（Forward-Time Centered-Space）
# ∂c/∂t + v ∂c/∂x = D ∂²c/∂x²
# +++++++++++++++++++++++++++++++++++++++++++++++

def reidai082():
    # ----- パラメータ設定 -----
    xmax  = 30.0        # 解析領域の長さ (m)
    dx    = 1.0         # 差分格子間隔 (m)
    dx2 = dx ** 2     # 差分格子間隔の二乗
    imax = int(xmax / dx) + 1    # 下流境界の格子点番号

    xo  = 10.0   # 濃度観測地点（m）
    ii1 = int(xo / dx)   # Python 0-based インデックス

    vel   = 0.1         # 実流速 (m/day)
    al    = 1.0         # 縦分散長 (m)
    tau   = 1.0         # 屈曲率 (-)
    dm    = 1.0e-5      # 分子拡散係数 (cm^2/s)
    # 分散係数 D (m2/day)
    # dm を m2/day に換算するために 1 cm2/s = 8.64 m2/day
    D     = al * abs(vel) + tau * dm * 8.64

    bc    = 1.0         # 上流境界の濃度 (mg/L)

    tend  = 300.0       # 計算時間 (day)
    dt    = 1.0         # 時間ステップ (day)
    nmax  = int(tend / dt)   # 計算時間のステップ数

    tjikan = 1.0        # 結果をファイルに出力する時間間隔（day）
    iout1 = int(tjikan / dt)   # 結果をファイルに出力する時間間隔のステップ数

    # ----- グリッド設定 -----
    x    = np.linspace(0, xmax, imax)

    # ----- 配列の用意 -----
    c  = np.zeros(imax)  # 新しい時刻の濃度
    co = np.zeros(imax)  # １ステップ前の濃度

    # ----- 係数の用意 -----
    r_diff = D * dt / dx2
    r_adv  = vel * dt / dx

    # 初期条件：上流境界だけ bc
    co[0] = bc

    # ----- 出力ファイル準備 -----
    out_filename = 'output.csv'
    with open(out_filename, mode='w', newline='') as f:
        writer = csv.writer(f)
        writer.writerow(['elapsed-time(days)', 'concentration(mg/L)'])
        # t=0 のデータ
        writer.writerow([f"{0.0:.1f}", f"{co[ii1]:.7f}"])

        # ----- 時間ループ -----
        time = 0.0
        iout  = 0
        for n in range(1, nmax+1):
            time += dt
            iout += 1

            # ---- 境界条件 ----
            co[0]      = bc      # 上流
            co[-1]     = 0.0     # 下流

            # ---- 内部格子の更新 (FTCS) ----
            # c[i] = co[i] + D*dt/dx^2*(co[i+1]-2*co[i]+co[i-1])
            #                  - vel*dt/dx*(co[i]-co[i-1])
            # スライスで一括計算
            # co[1:-1]：インデックス 1～N−2 の内部点
            # co[2:] ：インデックス 2～N−1
            # co[:-2] ：インデックス 0～N−3

            c[1:-1] = (co[1:-1]
                       + r_diff * (co[2:] - 2*co[1:-1] + co[:-2])
                       - r_adv  * (co[1:-1] - co[:-2]))

            # 更新
            co[:] = c[:]

            # ---- 出力 ----
            if iout == iout1:
                print(f"{time:6.1f} days")
                writer.writerow([f"{time:.1f}", f"{co[ii1]:.7f}"])
                iout = 0

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！お手軽に試すことができるようになりました。