waterPMLコマンドマニュアル

(The documentation of waterPML command)

Last Update: 2024/3/7

◆機能・用途(Purpose)

Maeda and Kumagai (2013)の差分法アルゴリズムを用いて 3次元空間内の弾性波動場を計算する。
Compute the elastic wavefield in a 3-D domain using the finite-difference algorithm of Maeda and Kumagai (2013).

このプログラムでは直交座標系$\posx=(x,y,z)=(x_0,x_1,x_2)$における線形弾性体の運動方程式 \[\begin{equation} \rho(\posx)\PartialDiff{V_i(\posx,t)}{t} =\sum_{k=0}^2 \PartialDiff{\tau_{ik}(\posx,t)}{x_k}+f_i(\posx,t) \hspace{1em} (i=0,1,2) \label{eq.motion} \end{equation}\] と構成則 \[\begin{equation} \PartialDiff{\tau_{ij}(\posx,t)}{t} =\sum_{p=0}^2 \sum_{q=0}^2 C_{ijpq}(\posx) \PartialDiff{V_p(\posx,t)}{x_q} \hspace{1em} (i,j=0,1,2) \label{eq.constitutive} \end{equation}\] を連立させて解く。ここで

$t$は時刻
$V_i(\posx,t)$は位置$\posx$、時刻$t$における速度の$x_i$方向成分
$\tau_{ij}(\posx,t)$は位置$\posx$、時刻$t$における応力の$(x_i,x_j)$方向成分
$f_i(\posx,t)$は位置$\posx$、時刻$t$において作用する力の$x_i$方向成分
$\rho(\posx)$は位置$\posx$における弾性体の密度
$C_{ijpq}(\posx)$は位置$\posx$における弾性定数

を表す。このプログラムでは等方弾性体 \[\begin{equation} C_{ijpq}(\posx) =\lambda(\posx)\delta_{ij}{\delta_pq} +\mu(\posx)(\delta_{ip}\delta_{jq}+\delta_{iq}\delta_{jp}) \label{eq.C} \end{equation}\] を仮定する。ここで$\lambda(\posx)$, $\mu(\posx)$は位置$\posx$におけるラメ定数を表し、 $\delta$はクロネッカーのデルタである。 $V_i(\posx,t)=0$, $\tau_{ij}(\posx,t)=0$を初期条件とし、ソースにおける等価体積力を$f_i(\posx,t)$で与えて $V_i(\posx,t)$, $\tau_{ij}(\posx,t)$の時間発展を計算する。解法の詳細は後述する。
This program solves the equation of motion (Eq. \ref{eq.motion}) and the constitutive law (Eq. \ref{eq.constitutive}) of a linear elastic medium in the cartesian coordinate system $\posx=(x,y,z)=(x_1,x_2,x_3)$, where

$t$ is time;
$V_i(\posx,t)$ is the $x_i$-component of the velocity at position $\posx$ and time $t$;
$\tau_{ij}(\posx,t)$ is the $(x_i,x_j)$-component of the stress at position $\posx$ and time $t$;
$f_i(\posx,t)$ is the $x_i$-component of the force exerted at position $\posx$ and time $t$;
$\rho(\posx)$ is the density of the elastic medium at position $\posx$; and
$C_{ijpq}(\posx)$ is the elastic constant of the medium at position $\posx$.

This program assumes an isotropic elastic medium (Eq. \ref{eq.C}), where $\lambda(\posx)$ and $\mu(\posx)$ are the Lame constants at position $\posx$ and $\delta$ is the Kronecker's delta. This program computes the temporal evolutions of $V_i(\posx,t)$ and $\tau_{ij}(\posx,t)$ with the initial conditions of $V_i(\posx,t)=0$ and $\tau_{ij}(\posx,t)=0$ and the equivalent body force of the source given by $f_i(\posx,t)$. Details of the computation method are described later.

◆ソースコード(Source code)

$YMAEDA_OPENTOOL_DIR/opentws/src/waterPML.c
$YMAEDA_OPENTOOL_DIR/opentws/src/waterPML_sub

◆使用方法(Usage)

コマンドライン引数でパラメータを指定する。パラメータの一覧を下表に示す。
Specify parameters by command-line arguments. The table below shows a list of parameters.

●「-」から始まらない引数 (Arguments not beginning with “-”)

このコマンドでは「-」から始まらない引数は存在しない。
This command does not have arguments not beginning with “-”.

●1つの「-」から始まる引数 (Arguments beginning with a single “-”)

このコマンドでは1つの「-」から始まる引数は存在しない。
This command does not have arguments beginning with a single “-”.

●「--パラメータ名=パラメータ値」の形式の引数 (Arguments of a form “--Parameter name=Parameter Value”)

「--パラメータ名=パラメータ値」の形式の引数は自由な順番で指定できる。「-」から始まらない引数の間に挿入しても良い。相反する指定がなされた場合には後の指定が優先される。
Arguments of a form “--Parameter name=Parameter Value” can be placed in an arbitrary order. They can even be inserted between arguments not beginning with “-”. In case of conflicting options being specified, the latter option has a higher priority.

(1)計算領域と格子セルに関するパラメータ (Parameters for the computational volume and grid cells)
waterPMLコマンドでは地球が球であることを無視して直交座標系で弾性波動場を計算する。 $x$軸を東方向、$y$軸を北方向、$z$軸を上方向 に取る。 $z$座標の原点($z=0$)は海水準とする。すなわち$z$座標は標高である。
The waterPML command computes the elastic wavefield with a cartesian coordinate system, ignoring the spherical geometry of the Earth. The $x$-, $y$-, and $z$-axes are taken east, north, and upward, respectively. The origin of the $z$-coordinate is taken at sea level. That is, a $z$-coordinate is equal to an altitude.

$x$, $y$, $z$軸に平行な直方体領域 (物理領域) を対象に計算を行う。これはメモリと計算時間の都合によるものであり、実際の媒質は物理領域の外まで続いているものと想定する。すなわち、物理領域の端は実際の媒質の端ではない。物理領域の端において現実には存在しない反射波が計算上で生じるのを避けるため、物理領域の外周に沿って波動を吸収する領域 (PML領域) を設ける。本マニュアルでは物理領域とPML領域を合わせた全領域を 計算領域 と呼ぶ(図1)。
The target region of the computation is a rectangular-solid volume (physical volume) parallel to $x$-, $y$-, and $z$-axes. This is because of limited memory and computation time; the actual medium is assumed to continue outside of the physical volume. That is, the end of the physical volume is not the end of the actual medium. To avoid waves reflected from the outer surfaces of the physical volume that are not present in the actual medium, the physical volume is surrounded by an artificial region that absorbs the waves (PML volume). Throughout this documentation, a computational volume means the total of physical and PML volumes (Fig. 1).

図1. 物理領域、PML領域、計算領域の関係。
Fig. 1. Relations of physical, PML, and computational volumes.

waterPMLコマンドでは物理領域を一様なサイズの直方体の小領域 (格子セル) に分割して計算を行う。網目状の線の交点としての「格子点」とは異なり、格子セルは直方体の領域を表すものである(図2)。
The physical volume is divided into small rectangular-solid sub-volumes (grid cells) of uniform size. A grid cell represents a region of rectangular-solid shape, different from a “grid node” that is an intersect of meshed lines (Fig. 2).

図2. (a)格子点と(b)格子セルの捉え方。このプログラムでは(b)の考え方を取る。
Fig. 2. Schematic illustrations of (a) a grid node and (b) a grid cell. This program uses the concept of (b).

物理領域の端は格子セルの表面となる。また、PML領域の厚さは格子セルのサイズの整数倍とする。すなわちPML領域の端も格子セルの表面である(図3)。
The end of the physical volume is on the surfaces of the grid cells. The thickness of the PML volume is set as an integer multiple of the grid cell size. Therefore the end of the PML volume is also on the surfaces of the grid cells (Fig. 3).

図3. 格子セルと物理領域、PML領域の関係。領域の塗りつぶし色は図1と同じ。
Fig. 3. Relation between grid cells and the physical and PML volumes. The fill colors for the volumes are consistent with those of Fig. 1.

以下のパラメータは物理領域、PML領域、格子セルを定義するためのものである。
The parameters below define the physical and PML volumes and grid cells.

パラメータ名 Parameter name	意味 Meaning	可能なパラメータ値 Allowed parameter values	デフォルト値 Default value
Nx	物理領域内の$x$方向の格子セルの数($N_x$)。 The number of grid cells along the $x$-direction in the physical volume ($N_x$).	3以上の整数。 An integer greater than or equal to 3.	省略不可 Cannot be omitted
Ny	物理領域内の$y$方向の格子セルの数($N_y$)。 The number of grid cells along the $y$-direction in the physical volume ($N_y$).	3以上の整数。 An integer greater than or equal to 3.	$N_x$
Nz	物理領域内の$z$方向の格子セルの数($N_z$)。 The number of grid cells along the $z$-direction in the physical volume ($N_z$).	3以上の整数。 An integer greater than or equal to 3.	$N_x$
Npmx	$x$軸の$-$側のPML領域の格子セルの数($N_{pml}^{-x}$)。 The number of grid cells in the PML volume on the negative side of the $x$-axis ($N_{pml}^{-x}$).	1以上の整数。 An integer greater than or equal to 1.	省略不可 Cannot be omitted
Npmy	$y$軸の$-$側のPML領域の格子セルの数($N_{pml}^{-y}$)。 The number of grid cells in the PML volume on the negative side of the $y$-axis ($N_{pml}^{-y}$).	1以上の整数。 An integer greater than or equal to 1.	$N_{pml}^{-x}$
Npmz	$z$軸の$-$側のPML領域の格子セルの数($N_{pml}^{-z}$)。 The number of grid cells in the PML volume on the negative side of the $z$-axis ($N_{pml}^{-z}$).	1以上の整数。 An integer greater than or equal to 1.	$N_{pml}^{-x}$
Nppx	$x$軸の$+$側のPML領域の格子セルの数($N_{pml}^{+x}$)。 The number of grid cells in the PML volume on the positive side of the $x$-axis ($N_{pml}^{+x}$).	1以上の整数。 An integer greater than or equal to 1.	$N_{pml}^{-x}$
Nppy	$y$軸の$+$側のPML領域の格子セルの数($N_{pml}^{+y}$)。 The number of grid cells in the PML volume on the positive side of the $y$-axis ($N_{pml}^{+y}$).	1以上の整数。 An integer greater than or equal to 1.	$N_{pml}^{-y}$
Nppz	$z$軸の$+$側のPML領域の格子セルの数($N_{pml}^{+z}$)。 The number of grid cells in the PML volume on the positive side of the $z$-axis ($N_{pml}^{+z}$).	1以上の整数。 An integer greater than or equal to 1.	$N_{pml}^{-z}$
x0	物理領域の$x$座標範囲の下端($x_0$)。[m] The lower limit of the $x$-coordinate range of the physical volume ($x_0$). [m]	任意の実数。 An arbitrary real number.	0.0
y0	物理領域の$y$座標範囲の下端($y_0$)。[m] The lower limit of the $y$-coordinate range of the physical volume ($y_0$). [m]	任意の実数。 An arbitrary real number.	0.0
z0	物理領域の$z$座標範囲の下端($z_0$)。[m] The lower limit of the $z$-coordinate range of the physical volume ($z_0$). [m]	任意の実数。 An arbitrary real number.	0.0
dx	格子セルの$x$軸方向のサイズ($\Delta x$)。[m] The size of grid cells along the $x$-direction ($\Delta x$). [m]	0.001以上の実数。 A real number greater than or equal to 0.001.	省略不可 Cannot be omitted
dy	格子セルの$y$軸方向のサイズ($\Delta y$)。[m] The size of grid cells along the $y$-direction ($\Delta y$). [m]	0.001以上の実数。 A real number greater than or equal to 0.001.	$\Delta x$
dz	格子セルの$z$軸方向のサイズ($\Delta z$)。[m] The size of grid cells along the $z$-direction ($\Delta z$). [m]	0.001以上の実数。 A real number greater than or equal to 0.001.	$\Delta x$

(2)時間に関するパラメータ (Parameters for time)
waterPMLコマンドでは時刻$t=t_{min}$からスタートして$t=t_{max}$まで、時間を順方向に$\Delta t$の刻みで動かしながら差分法により弾性波動場を計算する。以下はこれらの時間を定義するパラメータである。
The waterPML command computes the elastic wavefield in an ascending order of time at an increment $\Delta t$ from a start time $t=t_{min}$ to an end time $t=t_{max}$. The parameters below define these times.

パラメータ名 Parameter name	意味 Meaning	可能なパラメータ値 Allowed parameter values	デフォルト値 Default value
dt	計算の時間刻み($\Delta t$)。[s] The time stepping of the computation ($\Delta t$). [s]	$1.0\times 10^{-6}$以上の実数。 A real number greater than or equal to $1.0\times 10^{-6}$.	省略不可 Cannot be omitted
tmin	計算の先頭時刻($t_{min}$)。[s] The beginning time of the computation ($t_{min}$). [s]	$\Delta t$に非負の整数を掛けて得られる実数。 A real number that is a non-negative integer multiple of $\Delta t$.	0.0
tmax	計算の末尾時刻($t_{max}$)。[s] The end time of the computation ($t_{max}$). [s]	$\Delta t$に正の整数を掛けて得られる実数($>t_{min}$)。 A real number that is a positive integer multiple of $\Delta t$ and greater than $t_{min}$.	省略不可 Cannot be omitted

(3)PML領域での波の吸収に関するパラメータ (Parameters that define the absorption of the wave in the PML volume)
waterPMLコマンドではFesta and Nielsen (2003)に基づき、 PML領域内で微分方程式の$\PartialDiff{}{x_k}$ ($k=0,1,2$)を $\frac{i\omega}{i\omega+\alpha^k(|x_k-x_k^b|)}\PartialDiff{}{x_k}$ で置き換えることにより波を吸収する。ここで$\omega$は角周波数、$x_k^b$は最寄りの物理領域端の座標である。 $\alpha^k$は以下の式で定義される。
\[\begin{equation} \alpha^k(\xi)=A\frac{V_p(x_k^b)}{L_{pml}}\left(\frac{\xi}{L_{pml}}\right)^n \label{eq.alpha_k} \end{equation}\] ここで$V_p$はP波速度、$L_{pml}$はPML領域の厚さを表す。
In the waterPML command, the wave is absorved by replacing $\PartialDiff{}{x_k}$ ($k=0,1,2$) in the differential equations with $\frac{i\omega}{i\omega+\alpha^k(|x_k-x_k^b|)}\PartialDiff{}{x_k}$ in the PML volume following Festa and Nielsen (2003), where $\omega$ is the angular frequency and $x_k^b$ is the coordinate of the nearest boundary of the physical volume. Eq. (\ref{eq.alpha_k}) defines $\alpha^k$, where $V_p$ is the P-wave velocity and $L_{pml}$ is the thickness of the PML volume.

以下は(\ref{eq.alpha_k})式の$\alpha^k$を定義するパラメータである。
The parameters below define $\alpha_k$ in Eq. (\ref{eq.alpha_k}).

パラメータ名 Parameter name	意味 Meaning	可能なパラメータ値 Allowed parameter values	デフォルト値 Default value
PML_A	(\ref{eq.alpha_k})式の$A$の値。 The value of $A$ in Eq. (\ref{eq.alpha_k}).	0.0以上の実数。 A real number greater than or equal to 0.0.	10.0
PML_n	(\ref{eq.alpha_k})式の$n$の値。 The value of $n$ in Eq. (\ref{eq.alpha_k}).	0.0以上の実数。 A real number greater than or equal to 0.0.	2.0

(4)地形と水域に関するパラメータ (Parameters for the topography and water-filled regions)
waterPMLコマンドでは全領域を固体とする無限媒質での計算もできるが、地表面を与えてそこから上を真空とする、より現実的な問題設定での計算が可能である。
The waterPML command can compute the elastic wavefield in an infinite medium, where the entire region is assumed to be a solid medium. The command can also perform the computation in a more realistic situation, where the regions above a given ground surface is assumed to be the vacuum.

地表面は$(x,y)$座標ごと、あるいは緯度経度ごとの地表面の標高の値のリストとしてファイルで与える。このデータは数値標高モデルと呼ばれる。
The ground surface is given by a file that lists the altitude of the ground surface at each $(x,y)$-coordinate or at each latitude and longitude. This data is called a digital elevation model.

地形データは複数ファイルに分かれていても良いが、それらが全体として物理領域全域を包含している必要がある。地形データの定義域が物理領域の外にまではみ出している分には構わない。
The topography data can be separated into multiple files; however, the total of them must include the entire physical volume. There is no matter if the definition range of the topography data extends to the outside of the physical volume.

地形データはwaterPMLコマンド内部で格子セルに合わせて自動的に補間されるので、ユーザは格子セルに合わせたデータを用意する必要は無く、またそうすべきではない。ユーザがそのようなデータを用意した場合、 2重の補間が行われることで数値計算における地形表現の精度が落ちる場合がある。これを避けるため、元々の地形データをそのまま用いるべきである。
The topography data is automatically interpolated to adjust the grid cells within the waterPML command. Therefore the users do not need to, and should not, prepare the adjusted topography data. Use of such data may result in double interpolation, which may loose the accuracy of the topography expression in the numerical computation. To avoid it, the original topography data should directly be used.

waterPMLコマンドでは海や湖などの水域を計算に含めることができる。これには水面の標高や水の物性が書かれた設定ファイルを与える。海底・湖底の地形については地形データの一部として与えれば良い。地表面よりも水面の方が上にある領域が水域と見なされる。
The waterPML command can include water-filled regions (e.g., the sea and lakes) into the computation. This is realized by a configuration file in which the altitude of the water surface and physical properties of water are defined. Give the bathymetry of the sea and lake as a part of the topography data. Regions where the water surface is above the ground surface are regarded as water-filled regions.

waterPMLコマンドでは個々の格子セルに対して「固体」「水」「真空」のいずれかの属性を割り当てる。したがって地表面や水面は格子セルの境界面で近似されることになる(図4)。
The waterPML command assigns one of the “solid”, “water” or “vacuum” attributes to each grid cell. Therefore the ground and water surfaces are approximated at the boundaries of grid cells (Fig. 4).

図4. waterPMLコマンドにおける地形と水域の表現方法の模式図。
Fig. 4. Schematic illustration for the representations of the topography and water-filled regions in waterPML command.

以下は地形と水域に関するパラメータである。
Parameters for the topography and water-filled regions are shown below.

パラメータ名 Parameter name	意味 Meaning	可能なパラメータ値 Allowed parameter values	デフォルト値 Default value
topography_files	使用する地形データファイル名のリスト。 List of the file names of topography data used.	ファイル名をカンマ(,)で区切って書き並べた文字列。 A string composed of file names separated by commas (,).	省略時は地表面無し(無限媒質)での計算となる。 When this parameter is omitted, the computation is performed without the ground surface (i.e., an infinite medium).
topography_file_format	地形データファイルの書式。 The format of the topography data files.	latlon 緯度、経度、地表面の標高を組で与える。 Define the topography as combinations of a latitude, a longitude, and the altitude of the ground surface at that point. xy $x$座標、$y$座標、地表面の標高を組で与える。 Define the topography as combinations of an $x$-coordinate, a $y$-coordinate, and the altitude of the ground surface at that point. 3D ymaeda_opentoolsの3次元データ形式で与える。 Define the topography by a 3-D data of ymaeda_opentools. 詳細 (Detail)	3D
refN	緯度経度で与えられた標高データを直交座標に変換するための基準点(直交座標系の原点)の緯度。 The latitude of the reference point (the origin of a cartesian coordinate system) to convert the topography data, given in latitudes and longitudes, to the cartesian coordinates.	「度」単位の実数、または「度:分:秒」表記(度・分は整数)。 A real number in the degree unit, or a “degree:minute:second” representation in which the degree and minute are integers.	topography_file_format=latlonの場合は省略不可。その他の場合は使用されない。 Cannot be omitted if topography_file_format=latlon; in the other cases, this parameter is not used.
refE	緯度経度で与えられた標高データを直交座標に変換するための基準点(直交座標系の原点)の経度。 The longitude of the reference point (the origin of a cartesian coordinate system) to convert the topography data, given in latitudes and longitudes, to the cartesian coordinates.	「度」単位の実数、または「度:分:秒」表記(度・分は整数)。 A real number in the degree unit, or a “degree:minute:second” representation in which the degree and minute are integers.	topography_file_format=latlonの場合は省略不可。その他の場合は使用されない。 Cannot be omitted if topography_file_format=latlon; in the other cases, this parameter is not used.
topography_interpolation_method	地形データの補間に用いる方法。 The method of interpolation of the topography data.	3p 標高を計算したい地点の最寄りの3点のデータを用いて関数interpolate2D_3p により近似を行う。 Approximate the topography using function interpolate2D_3p with the data at the three points closest to the target point. mesh4p topography_file_format=3Dの場合にのみ選択可能である。標高を計算したい地点を囲む最小の長方形メッシュの角でのデータを用いて関数interpolate_3d_data_1p により近似を行う。 This choice is available only if topography_file_format=3D. Approximate the topography using function interpolate_3d_data_1p with the data at the four corners of the smallest rectangular mesh that bounds the target point.	topography_file_format=3Dの場合はmesh4p mesh4p if topography_file_format=3D その他の場合は3p 3p in the other cases
water_file	使用する水域の設定ファイル名。 The name of the configuration file for water-filled regions used.	ファイル名。 A file name. ファイルの書式 (File format)	省略時は水域無しでの計算となる。 When this parameter is omitted, the computation is performed without water-filled regions.

(5)地下構造に関するパラメータ (Parameters for the subsurface structure)
waterPMLコマンドでは任意の3次元P波速度・S波速度・密度構造を与えて計算を行うことができる。これらの地下構造情報は設定ファイルで与える。以下は地下構造の設定に関するパラメータである。
The waterPML command can perform the computation with arbitrary 3-D structures of the P-wave velocity, S-wave velocity, and density. These subsurface structural information are given by a configuration file. Parameters for the configuration of the subsurface structure are shown below.

パラメータ名 Parameter name	意味 Meaning	可能なパラメータ値 Allowed parameter values	デフォルト値 Default value
structure_file	使用する地下構造の設定ファイル名。 The name of the configuration file for the subsurface structure used.	ファイル名。 A file name.	省略不可 Cannot be omitted
structure_file_format	地下構造の設定ファイルの書式。 The format of the configuration file for the subsurface structure.	subdomain 物理領域を多面体の小領域に分割し、それぞれの小領域ごとにP波・S波速度と密度を指定する。 Divide the physical volume to polyhedron subdomains and specify the P- and S-wave velocities and density in each subdomain. 詳細 (Detail) layer 水平または地表面に平行な成層構造を与える。 Specify a stratified medium composed of layers that are horizontal or parallel to the ground surface. 詳細 (Detail)	subdomain

(6)地震波動ソースに関するパラメータ (Parameters for the seismic wave source)
waterPMLコマンドでは(\ref{eq.motion})式の$f_i(\posx,t)$をソース項として与えて波動場を計算する。以下は地震波動ソースに関するパラメータである。
The waterPML command gives the seismic wave source in the form of an equivalent body force $f_i(\posx,t)$ in Eq. (\ref{eq.motion}). Parameters for the seismic wave source are shown below.

パラメータ名 Parameter name	意味 Meaning	可能なパラメータ値 Allowed parameter values	デフォルト値 Default value
source_file	使用する地震波動ソースの設定ファイル名。 The name of the configuration file for the seismic wave source to be used.	ファイル名。 A file name. ファイルの書式 (File format)	省略不可 Cannot be omitted

(7)出力に関するパラメータ (Parameters that controll the outputs)
waterPMLコマンドで計算する量は位置$\posx$, 時刻$t$の関数としての速度場$V_i(\posx,t)$である。 waterPMLコマンドでは速度場を以下の2種類の方法で出力できる。

波形: 指定した位置$\posx$における、時間の関数としての速度場$V_i(t)$
スナップショット: 指定した時刻$t$における、位置の関数としての速度場$V_i(\posx)$

波形とスナップショットは両方出力することも片方だけを出力することも可能であるが、両方とも省略することはできない。なお、波形の出力においては厳密には観測点の位置ではなく観測点を含む格子セル中心点での速度場が出力される。
The waterPML command computes the velocity field $V_i(\posx,t)$ as a function of location $\posx$ and time $t$. The command can output the result in the following two ways:

waveforms, which are the velocity field $V_i(t)$ at specified locations as a function of time; and
snapshots, which are the velocity field $V_i(\posx)$ at specified time samples as a function of location.

The outputs can include the waveforms, the snapshots, or the both; at least one of them must be included in the outputs. The waveform outputs of the velocity field are not strictly at the given station locations but at the center of the grid cell that consists of each station.

waterPMLコマンドでは傾斜変動を出力することもできる。傾斜変動は鉛直速度成分$V_z$の水平方向への勾配として $T_x=\PartialDiff{V_z}{x}$, $T_y=\PartialDiff{V_z}{y}$ のように定義される。後述のパラメータtilt_outputは傾斜変動の出力方法を制御するパラメータであり、その意味を理解するためにはwaterPMLコマンド内部での格子セルを用いた速度場、傾斜変動、地形の表現方法を理解しておく必要がある。以下の3点がポイントである。
The waterPML command can also output the tilt changes that are defined as the horizontal gradiants of the vertical velocity $V_z$: $T_x=\PartialDiff{V_z}{x}$, $T_y=\PartialDiff{V_z}{y}$. A parameter tilt_output described later controlls the tilt outputs. To understand its meaning, users need to understand the internal representations of the velocity field, tilt, and topography on the basis of grid cells in the waterPML command. The following three points are essential.

waterPMLコマンドの内部では格子セルの上面と下面の中心点で $V_z$が定義されている。これを用いて格子セル中心点での傾斜変動は図5aのように計算される。
The waterPML command internally defines $V_z$ at the centers of the upper and lower surfaces of each grid cell. Using them, the tilt at the center of a grid cell is computed in a way schematically shown in Fig. 5a.
waterPMLコマンドでは個々の格子セルに「固体」「水」「真空」のいずれかの媒質を割り当てる。固体と真空のセルの境界面が地表面となる。すなわち、地形は格子セルの境界面によって階段状に近似される。
The waterPML command assigns one of the “solid”, “water”, or “vacuum” mediums to each cell. The boundary planes between the solid and vacuum cells constitute the ground surface. Therefore the internal expression of the topography is a stepwise approximation using the boundary planes of grid cells.
その結果、傾斜変動を計算したいセルの左右に固体の格子セルが存在しないケースが起こりうる(図4b-e)。このうち図4b, cのケースでは図の矢印で示した固体媒質内の$V_z$を用いて傾斜変動を定義できるが、図4d, eのケースでは傾斜変動を定義できない。
As a result, the target grid cell for the computation of tilt may not have adjacent grid cells to the right or left (Figs. 4b-e). In cases of Figs. 4b and c, the tilt can be computed from $V_z$ in the solid medium as illustrated by the arrows. In cases of Figs. 4d and e, however, the tilt cannot be defined.

図5. 傾斜変動の計算方法(青矢印)と計算できないケース。
Fig. 5. Computation method for tilt changes (blue arrows) and uncomputable cases.

このほか、計算で用いた内部データをチェック用に出力することもできる。以下は出力に関するパラメータである。
The outputs can also include internal data used in the computation for checks. Parameters that control the outputs are shown below.

パラメータ名 Parameter name	意味 Meaning	可能なパラメータ値 Allowed parameter values	デフォルト値 Default value
output_dir	計算結果の出力先ディレクトリ名。 The output directory name.	ディレクトリ名。既に存在するディレクトリ名を指定した場合には上書きされるので注意。 A directory name. Note that specifying an existing directory results in overwriting it.	.
波形の出力に関するパラメータ (Parameters for the output waveforms)
station_file	波形を出力する地点(観測点) のリストファイル名。 The name of a file that lists the locations where the waveforms are to be included in the outputs (stations).	ファイル名。 A file name. ファイルの書式 (File format)	省略時は波形を出力しない。 When this parameter is omitted, waveforms are not included in the outputs.
output_dt	出力する波形の時間刻み($\Delta t_{waveform}$)。このパラメータを指定することで計算の時間刻み$\Delta t$よりも粗いサンプルで出力することができる。但し単純に間引くだけであり、アンチエイリアスフィルターは適用されないことに留意。 The sampling interval of the output waveforms ($\Delta t_{waveform}$). By specifying this parameter, the time stepping the output waveforms can be more sparse than that of the computation $\Delta t$. However, note that the data samples are simply decimated without anti-alias filters.	$\Delta t$の自然数倍。 A natural number multiple of $\Delta t$.	$\Delta t$
スナップショットの出力に関するパラメータ (Parameters for the output snapshots)
snapshot_grid	スナップショットを出力する空間範囲。 The spatial extent of the output snapshots.	whole 計算領域全体でのスナップショットを出力する。 Output the snapshots in the entire computational volume. all 物理領域でのスナップショットを出力する。 Output the snapshots in the physical volume. specify:file=ファイル名 specify:file=file name スナップショットを出力する空間範囲と間隔をファイルで指定する。 Specify the spatial extent and interval of the output snapshots by a file. ファイルの書式 (File format) none スナップショットを出力しない。 Do not output the snapshots.	none
snapshot_place	スナップショットに含める量の定義位置。 waterPMLコマンドでは食い違い格子を用いるため、計算時の定義位置が変数毎に異なるが、それをそのまま出力するか、格子セル中心点での値に換算して出力するかをこのオプションで指定する。 The definition points of the quantities in the snapshots. The waterPML command a staggered grid, which means that the definition points are different among the variables. This option specifies whether to output them at their original definition points or at the center of grid cells.	c 計算結果を格子セル中心点での値に換算して出力する。 Output the results at the center of grid cells. d 計算時の定義位置で計算結果を出力する (但し、地表面・水面でのスナップショットは格子セル中心点での値となる)。 Output the results at the original definition points of the computation; however, the results at the center of grid cells are written in the snapshots on ground/water surfaces.	c
snapshot_dt	出力するスナップショットの時間刻み($\Delta t_{snapshot}$)。このパラメータを指定することで計算の時間刻み$\Delta t$よりも粗いサンプルで出力することができる。但し単純に間引くだけであり、アンチエイリアスフィルターは適用されないことに留意。 The sampling interval of the output snapshots ($\Delta t_{snapshot}$). By specifying this parameter, the time stepping the output snapshots can be more sparse than that of the computation $\Delta t$. However, note that the data samples are simply decimated without anti-alias filters.	$\Delta t$の自然数倍。 A natural number multiple of $\Delta t$.	$\Delta t$
その他の出力に関するパラメータ (Parameters for the other outputs)
output_parameters	計算で用いた格子セル毎の細かいパラメータの値をチェック用に出力するか否か。 Whether to output the detailed values of parameters for each grid cell used in the computation for check.	yes 出力する。処理時間や出力ファイルサイズが大きくなることに留意。 Output. This increases the processing time and output file size. no 出力しない。 Do not output.	no
output_stg	計算に用いた配列要素番号と空間内の位置の細かい対応関係をチェック用に出力するか否か。 Whether to output the detailed relations between the array indices and corresponding spatial locations used in the computation for check.	yes 出力する。出力ファイルサイズが劇的に大きくなることに留意。 Output. This significantly increases the output file size. no 出力しない。 Do not output.	no
tilt_output	観測点位置で傾斜変動を定義できない場合(図4d, e)の処理。 Treatment when the tilt cannot be defined at a station location (Figs. 4d and e).	strict プログラムをエラー終了する。 Finish the program as an error. none 計算は行うが、いずれの観測点においても傾斜変動は出力しない。 Continue the computation, but do not output the tilt at any stations. ignore 計算は行い、傾斜変動を定義できない観測点では傾斜変動を出力しない。 / Continue the computation, and do not output the tilt at stations where the tilt can not be defined. downward 計算は行い、傾斜変動を定義できない観測点ではその鉛直下方にある傾斜変動を定義できる格子セルでの傾斜変動を代わりに出力する。 Continue the computation, and for each station where the tilt change cannot be defined, output the tilt at a grid cell vertically downward where the tilt change can be defined.	strict

(8)計算の次元に関するパラメータ (Parameters for the computation dimension)
waterPMLコマンドでは3次元波動場を計算するが、 $y$方向に一様であることを仮定して $(x,z)$面内での2次元計算を行うこともできる。但しコーディングの都合上、本当の2次元計算ではなく3次元格子セルを用いた擬似的な2次元計算になっている。格子セルを3次元的に配置しても、 $y$方向に変化しない地形と地下構造を与え、速度の$y$方向成分$V_y$と応力成分$\tau_{xy}$, $\tau_{yz}$を初期値0から更新しなければ、他の速度場・応力場は$y$方向に一様な値となる。このような考え方で2次元計算が実装されているので、ユーザが意識する点として

$y$方向にも格子セルを配置する必要があること
地形と地下構造を$y$方向に一様な3次元構造として与える必要があること
ソースと観測点の$y$座標を揃える必要があること

に留意する。以下は計算の次元に関するパラメータである。
The waterPML command can not only compute 3-D wavefields and also 2-D wavefields assuming that any quantities do not change in the $y$-direction. However, this is not an actual 2-D computation but a pseudo 2-D computation in that the computation uses 3-D grid cells. The velocity and stress fields are uniform along the $y$-direction if the topography and subsurface structure are uniform along $y$ and the $y$-component of the velocity field $V_y$ and stress components $\tau_{xy}$ and $\tau_{yz}$ are not updated from their initial values of zero. Because the 2-D computation is implemented in this way, users must be aware of the following points:

users must assign grid cells along the $y$-direction;
the topography and 3-D structure must not change in the $y$-direction; and
the $y$-coordinates of the source and stations must be same.

Parameters that define the computation dimension are shown below.

パラメータ名 Parameter name	意味 Meaning	可能なパラメータ値 Allowed parameter values	デフォルト値 Default value
dimension	計算の次元。 The computation dimension.	2 2次元計算。 A 2-D computation. 3 3次元計算。 A 3-D computation.	3

(9)長時間の計算を補助するパラメータ (Parameters that support computation of long elapsed time)
スーパーコンピュータなど共用の高性能計算機では 1つのジョブに利用できる計算時間に制限がある場合がある。 waterPMLコマンドではそのような計算機を用いて 1度の計算では制限時間内に完了しない解析を複数回に分けて実行するための仕組みを用意している。
High-spec public computers (e.g., super computers) may have limits for the elapsed time of a single job. To adapt this situation, the waterPML command can divide an analysis into multiple jobs that should sequentially be executed.

その仕組みは以下の考え方に基づく。 waterPMLコマンドでは速度場と応力場の時間発展を解く。通常の計算では速度・応力の初期値を0として計算を行う。しかし計算が1度のジョブで完了しない場合、 1回目のジョブの最終ステップでの速度場・応力場全体をファイルに保存しておき、それらを初期値として2回目のジョブを走らせれば 1回目のジョブで行った計算の続きを行うことができる。同様にして3回目、4回目、…のジョブを走らせることができ、 1つのジョブでは制限時間内に完了しない長時間の波形を計算できる。
The idea is as follows. The waterPML command solves the temporal evolutions of velocity and stress fields. Usually, the initial values of the velocity and stress are set to 0. The division of an analysis into multiple jobs is realized by preserving the entire velocity and stress fields at the final step of the first job into a file, and restarting the computation in the second job using these preserved velocity and stress fields as the initial values. In the same way, the third, fourth, …, jobs can be executed. In this way, users can compute long waveforms than cannot be completed by a single job within the limited time.

以下はこのように計算を複数回のジョブに分割して行うためのパラメータである。
Parameters to realize these multiple jobs for an analysis are shown below.

パラメータ名 Parameter name	意味 Meaning	可能なパラメータ値 Allowed parameter values	デフォルト値 Default value
save_final	計算の最終ステップでの速度場・応力場をファイルに保存するか否か。 Whether to save the velocity and stress fields at the final time step into files.	yes 保存する。 Save. no 保存しない。 Do not save.	no
use_prev	速度場・応力場の初期値として前回の計算の最終ステップでの値を使用するか否か。 Whether to use the final values of the previous computation for the initial values of the velocity and stress fields.	yes 使用する。前回の計算において--save_final=yesオプションを用いている必要がある。 Use. This requires --save_final=yes option to have been used in the previous computation. no 使用しない。すなわち速度場・応力場の初期値を0とする。 Do not use; use 0 for the initial values of the velocity and stress fields.	no

◆動作(Behaviour)

(\ref{eq.motion})(\ref{eq.constitutive})式を連立させて弾性波動場の時空間発展を計算し、得られた速度場$V_i(\posx,t)$を波形またはスナップショットとして出力する。
Compute the spatiotemporal evolution of the elastic wavefield by simultaneously solving Eqs. (\ref{eq.motion}) and (\ref{eq.constitutive}), and output the resultant velocity wavefield $V_i(\posx,t)$ in the forms of waveforms or snapshots.

●出力されるファイル(Output files)

メインの計算結果は以下の2種類のファイルに出力される。パラメータによっては片方しか出力されない場合がある。
The main computation results are written into the following two types of files. Depending on parameters, only one of them may be created.
- 観測点での波形 (The waveforms at stations)
- スナップショット (Snapshots)
以下のファイルはパラメータsave_finalの値をyesにした場合に出力され、続きの計算を行う際に初期条件として用いられる。
The files below are created if the value of parameter save_final is yes; they are used to determine the initial condition for the continued computation.
- 最終時刻での速度場と応力場 (Velocity and stress fields at the final time step)
このほかに下記のファイルも出力することができる。これらのファイルは通常は不要であるが、プログラムのチェックや計算条件の技術的詳細を調べたい場合に使用する。
In addition, files listed below may be created. These files are usually not necessary, but may be used for checking the program and surveying technical details of computation condition.
- 実行したコマンドのログ (The log of the executed command)
- 使用したパラメータの値 (The values of parameters used)
- 異なる方法による位置の表現の変換表 (A conversion table for different expressions of locations)
- 半格子点の属性 (The attributes of half-grid nodes)
- 計算に使用した地表面の標高 (Ground surface elevation used in the computation)
- 計算に使用した水面の標高 (Water surface elevation used in the computation)
- 使用した地震波動ソースの情報 (Information on seismic wave sources used)
- 観測点の情報 (Information on stations)
- 観測点の周囲の地形 (Topography around stations)
- スナップショットの設定 (Configuration of snapshots)
- 密度の空間分布 (Spatial distribution of density)
- ラメ定数$\lambda$の空間分布 (Spatial distribution of Lame constant $\lambda$)
- ラメ定数$\mu$の空間分布 (Spatial distribution of Lame constant $\mu$)
- PML領域における吸収の強さ$\alpha^k$の空間分布 (Spatial distribution of the strength, $\alpha^k$, of absorption in the PML volume)
- 差分方程式の係数$c_0^k$の空間分布 (Spatial distribution of coefficient $c_0^k$ in difference equations)
- 差分方程式の係数$c_1^k$の空間分布 (Spatial distribution of coefficient $c_1^k$ in difference equations)
- 差分方程式の係数$c_2^{ijpl}$の空間分布 (Spatial distribution of coefficient $c_2^{ijpl}$ in difference equations)
- 等価体積力を位置の関数と時間の関数の積の重ね合わせで表したときの位置の関数の空間分布 (Spatial distribution of a position function when the equivalent body force is expressed as a summation of products of position and time functions)
- プログラムのデバッグのためのファイル (Files for debugging the program)

◆使用例と検証 (Examples and validations)

無限均質媒質中での計算 (Computation in a homogeneous infinite medium)
半無限媒質での計算 (Computation in a half space)
傾いた斜面での傾斜変動の計算 (Computation of tilt changes on a slope)
実際の3次元地形を用いた計算 (Computation using actual 3-D topography)
水域を用いた計算 (Computation using water-filled regions)

◆計算式とアルゴリズム (Formula and algorithm)

waterPMLコマンドの計算方法は以下の3種類のアルゴリズムを組合せたものである。
The computation method of waterPML command is the combination of the following three algorithms.

Ohminato and Chouet (1997)
任意の3次元地形と地下構造を取り入れた計算のための食い違い格子。
A staggered grid scheme to realize a computation including arbitrary 3-D topography and subsurface structure.
Festa and Nielsen (2003)
物理領域の端において反射波を吸収するためのPerfectly Matched Layer (PML)。
Perfectly matched layers (PML) to efficiently absorb the waves reflected from the boundaries of the physical volume.
岡本・竹中(2005)
Okamoto and Takenaka (2005)
海や湖などの水域を扱うアルゴリズム。
An algorithm to treat water-filled regions (e.g., sea and lakes).

このうちの1と2を組合せたものがMaeda et al. (2011)であり、更に3を組合せたものがMaeda and Kumagai (2013)である。 waterPMLコマンドで用いているのは速度-応力型の微分方程式であり、時間・空間方向ともに2次の差分を用いている。詳細を下記のリンクにて説明する。
Maeda et al. (2011) is the combination of algorithms 1 and 2, and Maeda and Kumagai (2013) is the combination of all three algorithms. The waterPML command uses differential equations in a velocity-stress form which are discretized using 2nd-order differences in both the temporal and spatial dimensions. See the links below for more detail.

微分方程式 (Differential equations)
離散化 (Discretization)
計算領域の外周での境界条件 (Boundary conditions on the outer surfaces of the computational volume)
地形・水域の表現方法と媒質境界での境界条件 (The representation methods for the topography and water-filled regions and the conditions on material boundaries)
地震波動ソースの与え方 (Defining the seismic wave source)

◆追加の情報 (Additional information)

ソースコードの概略(プログラマ向け) (Outline of the source code; for programmers)
ヘッダファイルと関数(プログラマ向け) (Header files and functions; for programmers)

◆引用文献 (References)

★ は本プログラムを用いた研究成果の発表時に必ず引用すべき論文を表す。
The star (★) indicates the references that must be cited whenever a research that used this program is published or presented.

Festa G, Nielsen S (2003) PML absorbing boundaries, Bull Seism Soc Am 93, 891-903. https://doi.org/10.1785/0120020098.
Gardner GHF, Gardner LW, Gregory AR (1974) Formation velocity and density — the diagnostic basics for stratigraphic traps. Geophysics 39(6), 770-780. https://doi.org/10.1190/1.1440465.
★ 水域を用いた計算の場合に引用のこと (Cite this article when your computation includes water-filled regions)
Maeda Y, Kumagai H (2013) Effects of water domains on seismic wavefields: A simulation case study at Taal volcano, Philippines, Earth Planets Space 65(2), 85-96. https://doi.org/10.5047/eps.2012.07.004
Maeda Y, Watanabe T (2023) Seismic structure and its implication on the hydrothermal system beneath Mt. Ontake, central Japan, Earth Planets Space 75, 115. https://doi.org/10.1186/s40623-023-01870-z
★ 水域を用いない計算の場合に引用のこと (Cite this article when your computation does not include water-filled regions)
Maeda Y, Takeo M, Ohminato T (2011) A waveform inversion including tilt: method and simple tests, Geophys J Int 184(2), 907-918. https://doi.org/10.1111/j.1365-246X.2010.04892.x
Maeda Y, Kato A, Terakawa T, Yamanaka Y, Horikawa S, Matsuhiro K, Okuda T (2015) Source mechanism of a VLP event immediately before the 2014 eruption of Mt. Ontake, Japan, Earth Planets Space 67, 187. https://doi.org/10.1186/s40623-015-0358-0
Mogi K (1958) Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them, Bull Earthq Res Inst 36(2), 99-134. https://doi.org/10.15083/0000033924
Ohminato T, Chouet BA (1997) A free-surface boundary condition for including 3D topography in the finite-difference method, Bull Seism Soc Am 87, 494-515. https://doi.org/10.1785/BSSA0870020494
岡元太郎, 竹中博士 (2005) 速度・応力型差分法での固体・流体境界の扱いについて, 地震第2輯, 57(3), 355-364.
Okamoto T, Takenaka H (2005) Fluid-solid boundary implementation in the velocity-stress finite-difference method, J Seism Soc Jpn 57(3), 355-364. (in Japanese with English abstract)
https://doi.org/10.4294/zisin1948.57.3_355
坂井孝行 (2012) 弾性体内部の変位とボアホール型傾斜計による傾斜変動量(第2報) 日本火山学会2012年度秋季大会, B3-01.
Sakai T (2012) Internal displacement of elastic body and tilt measured by borehole-type tiltmeter (2), Volcanol. Soc. Jpn. Fall Meeting 2012, B3-01. (in Japanese)
https://doi.org/10.18940/vsj.2012.0_82
武尾 (1985) 非弾性減衰を考慮した震源近傍での地震波合成 — 堆積層での非弾性減衰の効果について, 気象研究所研究報告, 36(4), 245-257.
Takeo M (1985) Near-field synthetic seismograms taking into account the effects of anelasticity — the effects of anelastic attenuation on seismograms caused by a sedimentary layer. (in Japanese with English abstract)
https://doi.org/10.2467/mripapers.36.245

waterPMLコマンド マニュアル