津波の数値シミュレーション

地学現象の多くは実験室ではほとんど再現できない，大きな時間・空間スケールのものが普通です。そのため，観測された事象の空間分布をモデル化し，時系列解析が行われますが，現在ではコンピューターをつかったシミュレーションが行われ，それによる理論と観測結果の比較研究や，予測による防災などにも利用され，非常に有力な手法となっています。ただし，中学高校の授業で取り上げるのは，やや高度な数学やコンピューターに関する知識が必要になるので難しいかもしれません。とはいえ，その原理など，ある程度理解しておくことも重要になっていると思います。昨年たまたま，教材化している例を，教員免許の更新講習で経験しましたので，十分理解できたとは言えないのですが，参考になればと思い紹介したいと思います。講習を受けたのは，2021年度のJpGU（日本地球惑星科学連合）主催による「数値シミュレーションで学ぶ津波の基礎」（講師；東京大学付属海洋教育センター准教授，丹羽淑博氏による）です。

これらの説明の内容を理解するには，微分方程式などに習熟していなければなりません。私も理解は十分とは言えませんが，要は一番下ような（黄色部分）の漸化式をコンピューターによって逐次計算させることができる（離散化させる）ということが，コンピューターシミュレーションのポイントであると思います。講習では，プログラミング言語として，Fortran（gfortran)が使われ，アニメーションにはGnuplotというソフトが使われましたが，私の使い慣れているPython に実習で配布されたFortranのコードを，移植してアニメーションにしたものを順に載せておきます。

津波の伝わる速さは水深ｈの平方根に比例します。沖合に戻る波の速さの変化に注目

以上の参考資料：
丹羽淑博・佐藤俊一（2021)：数値シミュレーションで学ぶ 津波の物理の基礎．
日本地球惑星科学連合2021年度教員免許状更新講習「数値シミュレーションで学ぶ津波の基礎」（2021年12月27日実施）
https://www.cole.p.u-tokyo.ac.jp/common/img/mt/events/pdf/20110827_niwa.pdfも参考にしました。

　以下は，講習で配布されたgfortranのコードをPythonで書き直したものです。

Pythonによる計算プログラム（Tsunami01）
　差分近似式に対応する変数
過去        U(t—△t)   →    ub
現在        　U(t)　 　→　   up
未来        U(t＋△t)   →   　uf
過去        Z(t—△t)　 →     zb
現在        　Z(t)　     →　  zp
過去        Z(t＋△t)    →     zf

#-------------------------------------------------------------------------------
# Name:        Tsunami01
#
# Purpose:     水槽内の水の水位と流速の変化（津波のシミュレーション）
#
#-------------------------------------------------------------------------------
import math
import numpy as np
import matplotlib.pyplot as plt
# from matplotlib.animation import FuncAnimation
from matplotlib.animation import ArtistAnimation
#---グリッド分割数-------
iq = 500
#----パラメータ値を設定-------------
up = [0]*(iq+1)
zp = [0]*(iq+1)
uf = [0]*(iq+1)
ub = [0]*(iq+1)
zb = [0]*(iq+1)
zf = [0]*(iq+1)
xlist = [[0 for i in range(iq+1)] for j in range(240)]
zlist = [[0 for i in range(iq+1)] for j in range(240)]
ulist = [[0 for i in range(iq+1)] for j in range(240)]
depth = 200.           #水深 (m)
width = 100.*1000.     #モデル海洋の幅 (m)
grav = 9.8             #重力加速度  (m/s^2)
#------------------------------
time_end = 2.*60.*60.  #計算時間(sec)
time_out = 30.         #データ出力時間間隔(sec)
#----------------------------
dx = width/iq      #空間グリッドサイズ
dt = 1             #時間ステップサイズ
#--初期値の設定-------
for i in range(iq+1):
    zp[i] = math.exp(-(i-iq/2)**2/(iq/20)**2)

nend = int(time_end/dt)  #時間ステップ数
nout = int(time_out/dt)  #データ出力ステップ間隔
nc = 0
#----------時間ループ n=0,nendまで繰り返す----
for n in range(0,nend):

#**************************************************************
#-----最初のステップのみ現在ステップから未来ステップ値を計算（前方差分）----
    if n==0 :
        for i in range(1,iq):
            uf[i] = up[i] - grav*(dt/dx)*(zp[i+1] - zp[i])

        for i in range(2,iq):
            zf[i] = zp[i] - depth*(dt/dx)*(up[i] - up[i-1])

#-----過去ステップ値と現在ステップ値から未来ステップ値を計算（中央差分）----
    if n>=1 :
        for i in range(1,iq):
            uf[i] = ub[i] - 2*grav*(dt/dx)*(zp[i+1] - zp[i])

        for i in range(2,iq):
            zf[i] = zb[i] - 2*depth*(dt/dx)*(up[i] - up[i-1])

#---境界条件----------
    uf[1] = 0.
    uf[iq-1] = 0.
#----計算の安定化のためのおまじない（Asselin filter)---
    if n>=1 :
        eps=0.01
        for i in range(iq):
            up[i] = up[i] + eps*(uf[i] - 2*up[i] + ub[i])

        for i in range(iq+1):
            zp[i] = zp[i] + eps*(zf[i] - 2*zp[i] + zb[i])

 #---noutステップごとに図のデータとして保持 --------------
    if n % nout == 0:
        for i in range(iq+1):
            xlist[nc][i] = dx*(i-2)/1000
            zlist[nc][i] = zp[i]
            ulist[nc][i] = up[i]

        nc += 1

#------ステップを進める。現在ステップ値を過去ステップ値に、
#      未来ステップ値を現在ステップ値に----
    for i in range(iq):
        ub[i] = up[i]
    for i in range(iq+1):
        zb[i] = zp[i]
    for i in range(iq):
        up[i] = uf[i]
    for i in range(iq+1):
        zp[i] = zf[i]

#---------アニメーション
fig = plt.figure()
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
plt.subplots_adjust(hspace=0.4)
artists = []
for n in range(nc):
    zline, = ax1.plot(xlist[n], zlist[n],"blue")
    ax1.set_title("surface elevation (m)  ")
    vline, = ax2.plot(xlist[n], ulist[n],"blue")
    ax2.set_title("horizontal velocity (m/s)")

    artists.append([zline,vline])

anim = ArtistAnimation(fig, artists, interval = 100, repeat_delay = 600)

# anim.save("Tsunami1.gif", writer="imagemagick")
plt.show()