Rogiers et al.(2022), Geophys. J. Int.；降下テフラの噴出量推定に関するメモ書き
研究の関係上降下テフラの噴出量推定に関する論文は少し漁っているつもりですが、今春こんな論文が出ていたので今回はこれを叩き台に幾つか感想という名の戯言を書き連ねていきたいと思います。
Rogiers et al.(2022), Geophys. J. Int.
Estimating tephra fall volume from point-referenced thickness measurements
降下テフラの噴出量推定手法は、今は廃れたように思える結晶法の系譜（日本の火山学界隈が何故か未だによく使う早川法もここに含めることとする）を除くと
①等層厚線の層厚-(面積)^(1/2)図を描いて経験式でフィッティングする
②移流拡散モデルに基づくインバージョン
の２つに大別されると思いますが、②は給源強度関数や風速場に敏感だったり、パラメータが多かったり、粒度分析をしにくい保存状態が悪いテフラ層には不向きだったり…と今僕がやっている有史以前の噴火に関しては扱いにくい面がありそうなので、今のところは①に着目しています（と言いつつ実は単にC言語に触れたくないだけ…）。噴出率とか他のパラメータとも関連付けられてそれはそれで楽しそうだし、いずれはTephra2なりWTも弄ってみたいですね。
①の手法に関して頻繁に使われる有名どころは以下のようなものでしょうか。
Pyle(1989), Bull. Volcanol.
The thickness, volume and grainsize of tephra fall deposits
https://doi.org/10.1007/BF01086757
Fierstein and Nathenson(1992), Bull. Volcanol.
Another look at the calculation of fallout tephra volumes
https://doi.org/10.1007/BF00278005
Bonadonna and Houghton(2005), Bull. Volcanol.
Total grain-size distribution and volume of tephra-fall deposits
https://doi.org/10.1007/s00445-004-0386-2
Bonadonna and Costa(2012), Geology.
Estimating the volume of tephra deposits: A new simple strategy
https://doi.org/10.1130/G32769.1
Pyle(1989)とFierstein and Nathenson(1992)は指数関数、Bonadonna and Houghton(2005)はべき関数、Bonadonna and Costa(2012)はWeibull関数をそれぞれ使っています。最近は4番目のものを使った論文が多く出始めている印象。
この4手法、噴出量推定値に結構な差が出ます。下は卒論やJpGU2022で僕が某テフラに適用した結果です。
さて、ここまで色々話してきましたが、この話をとあるゼミでしたところある先輩からこんな指摘が飛び出してきました。
「その推定値の信頼区間はどれくらいなんですか？そもそもその等層厚線に信頼性はどれくらいあるんですか？」
…正直言葉に詰まりました。噴出量もそうですが、火山地質学の世界では数値に信頼区間を付ける習慣がほとんど無いように見えます。火山学の外の人から言われるまでそんなに気にしていませんでしたが、確かにまずいですね。
というわけで、今年度初めくらいは火山地質学における噴火パラメータ推定値に信頼区間を付ける論文をいくらか読んでいました。
Biass et al.(2014), Stat. Volcanol.
TError: towards a better quantification of the uncertainty propagated during the characterization of tephra deposits
https://doi.org/10.5038/2163-338x.1.2
Biass et al.(2019), J. Appl. Volcanol.
A step-by-step evaluation of empirical methods to quantify eruption source parameters from tephra-fall deposits
https://doi.org/10.1186/s13617-018-0081-1
どちらもBonadonna一派のところですね。２つありますが、手法の要諦は同じです。要は、今までと同じように手書きで等層厚線を描いて、その層厚や囲まれる面積の信頼度を適当に決めてやれば、後はフィッテング時の誤差伝播を考えれば噴出量推定値にも信頼区間らしきものが出てきます。余談ですがMatlabコードが配布されているので僕みたいなプログラミング弱者にもとっつきやすい。
しかし、僕は前々から思っていることがありました。
「そもそも等層厚線を人間が描いている以上そこに何らかの恣意性が入り込んでしまうのでは？」
そんな一例として良いのは下の論文でしょうか。阿蘇4の噴出量推定に関するものですが、本論文中で示されている阿蘇4火山灰の等層厚線に関してはこの疑問点が常に付きまとっているように見えます。
Takarada and Hoshiizumi(2020), Front. Earth Sci.
Distribution and Eruptive Volume of Aso-4 Pyroclastic Density Current and Tephra Fall Deposits, Japan: A M8 Super-Eruption
https://doi.org/10.3389/feart.2020.00170
等層厚線に関してはEngwellらがCubic-Bスプライン曲線を用いて色々やっていたのをちらっと見ましたが、コードが配られていなかったのでそのまま読まなかったことにしています。プログラミング弱者はつらい。
そんな折、指導教員からある論文を紹介されました。それが冒頭に述べたRogiers et al.(2022)という訳です。
Rogiers at al.(2022)の手法の要諦は以下のような感じです。とは言っても、等層厚線描画時のアルゴリズムの仔細に関してはほとんど理解出来ませんでしたが…
層厚地点データとその信頼度の付与
→Box-Cox変換による層厚データの正規分布化
→ Delaunay triangulationを用いた仮想地点の創出
→楕円をbasic functionとしてLASSOを用いたフィッティングによる等層厚線の描画（この箇所は特に難解で最も理解できなかった）
→積分して噴出量を算出
かなり難解なシロモノですが、等層厚線描画時の恣意性を排除していること、層厚データそのものを確率変数化して厳密に統計学的議論を試みているのは魅力的です。
さて、本論文ではR言語によるコードが配られています。というかこれが僕の興味を惹いた主因。今回はこのコードで鬼界アカホヤ火山灰の噴出量を推定してみます。これ研究室ゼミで発表はしたんですが指導教員様の反応が芳しくなかったので此処で供養。
Supplementary Dataのところにzipファイルが置かれているので解凍して中身を見てみるとこんな具合。
pdfファイルは筆者らが実際に阿蘇4テフラに適用した結果の図なのでこれを使う上では要りません。削除。必要なのはtephra4_3, randomBasisと層厚データ（論文中では20210805_Aso4_Tephra）です。
層厚データは一定のフォーマットで書かれないと上手く読み取ってくれません。今回は鬼界アカホヤ火山灰の噴出量推定ということで、僕が今までやったフィールド調査や諸々の論文や遺跡発掘調査報告書から取り敢えず層厚データを引っ張ってきて20210805_Aso4_Tephraを書き換えていきます。
Identification_uncertainyやThickness_uncertaintyが信頼度付与の部分です。
次にコードを書き換えます。弄るのはtephra4_3の方です。
元のコードは阿蘇４なので、取り敢えずこれを全部鬼界アカホヤに書き換える。
給源火口の位置を打ち込みます。
pseudo-thicknessという手描き層厚線データを打ち込んでやります。これはアルゴリズムを走らせる際に実際の層厚データが無い地域を変な風に扱わないようにするためのものです。
配られているコードでは常に同じ図を吐くように乱数を調整していますが、実際に使う上ではエラーの元なのでコメントアウト（44行目）。これ本当にやっていいのかは正直分かりませんが…
これでコードの方の準備は多分大丈夫です。
次にこのコードを走らせます。次のパッケージをインストールする必要があります。
data.table
magrittr
deldir
polyclip
gmlnet
sp
他にある場合はコンソール画面で教えてくれます（適当）。
では動かしていきましょう。コマンドは
source("tephra4_3.R")
です。
止まるまでに結構時間がかかりますが、最終的にこんな画面が出てきます。
見かけ噴出量の95％信頼区間は98-140㎦と推定されました。この値は町田・新井(2003)による推定値（>100㎦)と大体同じオーダーですね。描かれた等層厚線図は下の通りです。
一般的に等層厚線図は層厚データの最大値を囲むように描かれることが多いですが、見ても分かるようにRogiersのアルゴリズムに基づいて描いたそれはそのようになっていません。そのため、上で示した95％信頼区間は、これまでのように人間のある意味「勘」によって描かれた等層厚線図に基づく推定値よりも小さいものとなっていることが想定されます。
さて、ここまでRogiersのアルゴリズムで遊んできましたが、では今までのように人間が描いた等層厚線をフィッティングする手法と比較した際に、この手法が客観性はともかく正確性の意味で本当に有利なのか？というのは僕には正直分かりません。理由としては、層厚データが原層厚を中心とした正規分布に従う根拠は無く、削剥によって薄くなることは頻繁にあっても、著しく厚くなることはそれほど多くないと考えられるからです。
結局、噴出量推定は手法やフィッティング関数によって数倍程度は変わってきてしまうのが普通なので、Rogiersらも本論文の中でしつこく言っているように噴出量を１つの値として示すのではなく、（出来れば複数の手法で推定しつつ）ある程度の幅で示すのが望ましいのでしょう（適当）。
以上、最後は投げやりですがRogiers+2022のコードで遊んでみた結果でした。
