火山灰追跡モデルPUFFをPythonで実装してみた
1. 飛行機にとってクリティカルな火山灰
火山の噴火は、空を飛ぶ航空機にとってクリティカルな現象の1つです。噴煙に入ってしまうと視界が悪くなるばかりでなく、火山灰によってコックピットの窓ガラスに傷が付くなど、機体の損傷を引き起こしてしまいます。
特に問題なのは、エンジンに入ってしまうと熱で溶けた火山灰がエンジンに付着し、最悪の場合はエンジン停止に陥ってしまうことがあります。幸いにして火山灰が原因で墜落したケースはありませんが、一時的に全エンジンが停止して緊急着陸を行った（エンジンは着陸前に再始動した）ケースは過去に発生しています。
このように火山灰は航空機の安全運航に大きな影響があるため、火山灰の拡散予測は気象庁の重要な業務の1つとなっています。
（参考サイト：東京航空路火山灰情報センター｜Tokyo VAAC）
さて、タイトルにもある火山灰追跡モデルPUFFは、火山灰の拡散予測モデルの1つです。筑波大学の田中博教授がUAF (University of Alaska Fairbanks) にいらっしゃったころ、1990年代前半に考案されました。PUFFは現在でもアンカレッジVAACなどで現役で活躍しており、実績あるモデルです。
（参考サイト：PUFF model）
今回ちょっとした縁があり、ある火山の噴火の航空路への影響を調査するため、PUFFを自分で実装することになりました。そこで自分自身の整理もかねて、PUFFについてnoteに記述したいと思います。
なお記述する内容のほとんどが、「6. 参考文献」に記載した論文をもとにしていますので、詳細を知りたい方は論文の方も参照して下さい。
2. 火山灰追跡モデルPUFF
最初にシミュレーション例をご紹介します。事例についての説明は「5. シミュレーション結果（例）」を参照して下さい。
PUFFのユニークな点として「ラグランジュ形式で火山灰の動きをシミュレーションする」という特徴があります。つまり数値予報のように各格子点上の火山灰濃度を時間積分して求めるのではなく、火山灰粒子の動きを1個1個それぞれ計算するという手法を取ります。
この手法だと、数値予報が計算した気象場を利用して事後的に火山灰の動きだけをシミュレーションできるため、噴火が発生したら比較的短時間で拡散予測を行ってエアラインへ情報提供できる、というメリットがあると私は理解しています。
PUFFの構造を簡単に書くと、次のようになります。
火山灰の動き ＝ 風による輸送 ＋ 乱流による拡散 ＋ 重力落下
風による輸送はその名の通り、風に流されて輸送される効果を表す項です。火山灰の未来の動きを予測したい場合はGSMなどの数値予報モデルの風予測値を使い、過去の噴火の火山灰拡散シミュレーションをしたい場合はJRA55などの再解析データを使うことで、風による輸送を計算できます。
ここで鉛直方向については、風の鉛直成分というデータはなく、鉛直P速度があっても重力落下に比べて小さいことから、鉛直方向の風の輸送は考えない仕様になっています。
またGSMでもJRA55でも、GPVデータの時間・空間解像度が粗いという問題があります。PUFFのオリジナルの仕様では時間ステップを5分で取りますが、GSM全球域は6時間間隔・GSM日本域は3時間間隔（※気圧面データ）・JRA55は6時間間隔となっており、GPVの時間間隔がかなり荒いです。そこでCubic Splineで5分間隔に内挿する、という方法を取ります。
空間についても同様に内挿します。私の実装の場合、鉛直方向についてはCubic Splineで100m間隔に内挿します。水平面についても、火山灰の位置に応じてCubic Splineで内挿した風予測値を取得します。
なお個人的には、高い計算コストを掛けてCubic Splineで内挿して、それに見合った効果があるのかな？とも思います。線形内挿でも良いように思いますが、一応オリジナルの仕様にそってCubic Splineで内挿するように実装しています。
乱流による拡散は、GSMやJRA55などのGPVでは表現できない気流の乱れによって、火山灰が拡散される効果を表す項です。正規分布に従って発生させた乱数を拡散速度として使い、風と同様に火山灰の移動量を計算します。拡散係数等のパラメータについては参考文献[1]などの論文を参照しています。
重力落下もその名の通り、重力によって火山灰が落下する効果を表します。落下速度は改良型ストークスの法則で求めます。具体的な式は、参考文献[1]などの論文を参照しています。
なお鉛直風をゼロと扱うので、鉛直方向の動きは乱流による拡散と重力落下によって決まります。私がざっと計算させてみたところ、大雑把にいって火山灰の粒径が10μm程度かそれ以下になると、重力落下を打ち消して乱流拡散のみで高度を上げることがあります。
この他、火山灰の拡散シミュレーションを本格的にやる場合は、雨で洗い流される効果や、一度地面に降灰した火山灰が風で舞い上がる効果なども考えるようですが、私の実装ではそこまでの効果は取り入れていません。
3. 噴煙の初期設定
火山灰の動きの計算は上記の通りですが、この他に噴煙についての初期パラメータが色々あります。
まず噴煙トップの高度。通常は衛星画像解析やライブカメラなどの観測を用いて噴煙トップの高度を推定しますが、不明な場合や事前に噴火の影響を見積もるためにシミュレーションする場合は、「大は小を兼ねる」的に高めの高度に設定して計算させるようです。
私の実装では、噴火開始時の噴煙トップの高度をパラメータとして設定し、時間ステップごとに一定の割合で逓減させるような仕様にしています。
次に火山灰粒子の数。当たり前ですが、実際の噴火で噴出される全ての火山灰粒子の動きを計算することは非現実的なので、初期の火山粒子数や同時に計算する最大粒子数を限定してシミュレーションします。そして1つの火山灰粒子で、ある程度の質量の火山灰の動きを代表させて計算するようです。
私の実装では噴火開始時の火山灰粒子数をパラメータとして設定し、時間ステップごとに一定の割合で逓減させるような仕様にしています。
火山灰粒子の粒径の設定も必要です。対数正規分布で乱数を発生させて粒径を与えます。ここで対数正規分布の平均と標準偏差の設定は、各論文を見ると噴火事例ごとに異なるようです。詳しい値は参考文献を参照して下さい。
そして噴出した火山灰の初期の高度分布。一番簡単なのは火口から噴煙トップまで一様に分布させる方法で、実際の噴煙の状態が不明な場合や、事前に噴火のシミュレーションをする場合などに使われるようです。
ただ実際の噴火では、噴煙トップに多くの火山灰が集まりますので、そのような高度分布になるように乱数を発生させる手法もあります。
さらに、火口から噴出した火山灰は大きな鉛直速度を持って上昇するため、噴煙トップで水平方向に拡散する効果も組み込むことがあるようです。具体的には、噴煙の高度と火山からの水平距離に応じて、水平方向の拡散係数を調整します。例えば高度が高く、火山に近い位置にある火山灰は、乱流による拡散よりも水平方向に大きな拡散速度を持たせます。詳細は参考論文[1]などを参照して下さい。
4. プログラミング言語
私自身が最近使い慣れていて、GPVデータの処理も可能なPythonで実装しています。計算速度を重視するならC系の言語やFortranの方が良いのかもしれませんが、本業プログラマーではない私でも扱いやすいPythonでやってみました。
PythonでPUFFを実装するって、私くらいでしょうか（笑）
処理速度は、どうもCubic Splineによる内挿に時間が掛かってそうですが、シミュレーション後の可視化も含めて、個人持ちのMacBook Proでも十分な速度で計算できている印象です。
いずれソースコードは公開するかもしれませんし、しないかもしれません。とりあえず今はまだ開発中でもあるし非公開です。
5. シミュレーション結果（例）
PUFFによるシミュレーション結果は、5分間隔で火山灰粒子一つ一つを地図上にプロットした画像を作成し、さらにそれらをつなげた動画まで作るようにしています。以下にその一例の動画を紹介します。
この動画は参考文献[1]と同じく、2014年2月13日(UTC)に発生したインドネシアのKelud Volcanoの噴火について、PUFFで火山灰拡散シミュレーションを行ったものです。
参考文献[1]とは使っている気象データが異なるなどの理由により、全く同じシミュレーション結果にはなっていませんが、論文に掲載された図とだいたい同じような計算ができていると思います。
6. 参考文献
[1] H. L. Tanaka, M. Iguchi, and S. Nakata 2016: Numerical simulations of volcanic ash plume dispersal from Kelud volcano in Indonesia on 13 February 2014. J. Disaster Res. 11, 1, 31-42.
[2] H. L. Tanaka and M. Iguchi 2016: Numerical Simulation of Volcanic Ash Plume Dispersal from Kuchinoerabujima. J. Natural Disaster Sci., 37, 2,79-90.
[3] H. L. Tanaka and M. Iguchi 2019: Numerical Simulations of Volcanic Ash Plume Dispersal for Sakura-jima using Real-time Emission Rate Estimation, Journal of Disaster Research, 14, No.1, 160-172, 2019
[4] 田中博 2003: リアルタイム火山灰追跡モデル(PUFF)を用いた予報実験
[5] 邉見萌乃 2011: 火山灰追跡モデルPuffを用いた2010年4月のアイスランド火山灰の輸送拡散数値実験
さいごに
せっかくソフトウェアを作ったので、何かビジネスに活用できないかなと考えています。もしPUFFを使った火山灰拡散シミュレーションや、過去の噴火の再現などに興味のある方、下記よりお気軽に問い合わせいただけましたら幸いです。
最後まで読んでいただき、ありがとうございました。
この記事が気に入ったらサポートをしてみませんか？