waterPMLコマンドで用いている計算式

(4) 地形・水域の表現方法と媒質境界での境界条件

(Formula used in waterPML command; (4) The representation methods for the topography and water-filled regions and the conditions on material boundaries)

waterPMLコマンドでは個々の格子セルに対して 「固体」「水」「真空」のいずれかの属性を割り当てる。 したがって地表面や水面は格子セルの境界面で近似されることになる(図1)。 以下では「固体」「水」「真空」の属性別の取り扱い、 およびそれらの境界の取り扱いについて述べる。
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. 1). Below, treatments for each attributes (“solid”, “water” and “vacuum” cells) and their boundaries.


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

4-1. 格子セルの中心点での物性値 (Physical properties at the centers of grid cells)

4-1-1. 固体のセル (Solid cells)

地下構造の設定ファイル (小領域分割 または 層構造) において地下の全ての地点の 密度\(\rho\)、P波速度\(V_p\)、S波速度\(V_s\)が定義される。 これらを元に物理領域内の固体の格子セル中心点での密度とラメ定数 \(\lambda=\rho(V_p^2-2V_s^2)\), \(\mu=\rho V_s^2\) を設定する。
A configuration file for the subsurface structure (either subdomain or layer) defines the density \(\rho\), P-wave velocity \(V_p\), and S-wave velocity \(V_s\) over the entire subsurface. Based on these configurations, the density and Lame constants \(\lambda=\rho(V_p^2-2V_s^2)\) and \(\mu=\rho V_s^2\) at the center of each grid cell of solid type in the physical volume are determined.

4-1-2. 水のセル (Water cells)

水域の設定ファイル において水の密度\(\rho_w\)と音速\(a_w\)が定義される。 音速はP波速度に対応し、S波は水中を伝播しない。 そこで、物理領域内の水の格子セル中心点での密度を\(\rho_w\)、 ラメ定数を\(\lambda=\rho_w a_w^2\), \(\mu=0\) と設定する。
A configuration file for water-filled regions defines the density \(\rho_w\) and sound velocity \(a_w\) of water. The sound velocity corresponds to a P-wave velocity, while the S-wave does not propagate in the water. Therefore, the density \(\rho_w\) and Lame constants \(\lambda=\rho_w a_w^2\) and \(\mu=0\) are set at the center of each grid cell of water type in the physical volume.

4-1-3. 真空のセル (Vacuum cells)

真空での計算は必要ない。 そこで、真空の格子セルにおいては密度やラメ定数を設定せず、 真空であることを表すフラグのみを設定する。
Computation is not needed in the vacuum. Therefore, densities and Lame constants are not defined in vacuum cells. Instead, a flag is defined to indicate that these cells are vacuum.

4-2. 格子セルの境界での物性値 (Physical properties on the boundaries of grid cells)

での密度とラメ定数が計算に必要になる。 周囲の格子セル中心点(真空やPML領域を除く)での密度やラメ定数を平均した値を これらの地点での密度、ラメ定数として用いる。
Densities and Lame constants are required at:
They are computed by averaging the densities and Lame constants at the centers of the surrounding cells (excluding vacuum cells and cells in the PML volume).
なお 固体のセルと水のセルの境界(水底)においては ラメ定数\(\mu=0\)とする(平均ではない)。 この点を除けば水のセルも固体のセルと同様に扱う。 この方法が水域の取り扱いとして適切であることが 岡本・竹中(2005)で示されている。
Water cells are treated in the same manner as solid cells except that the Lame constant \(\mu\) is set to zero (not an average) on the boundary between solid and water cells (water floor). Okamoto and Takenaka (2005) showed the validity of this treatment.

4-3. PML領域内での属性と物性値 (Attributes and physical properties in the PML volume)

PML領域内での格子セルの属性(固体/水/真空の区別)や密度とラメ定数の値は 最寄りの物理領域内の地点に合わせる。 The attribute of each grid cell (the distinction of solid, water, or vacuum cells) and the density and Lame constants in the PML volume are set to be equal to those at the nearest location in the physical volume.

4-4. 地表面や水面の境界条件 (Boundary conditions on the ground and water surfaces)

地表面(固体セルと真空セルの境界)や水面(水セルと真空セルの境界)においては 自由表面境界条件(トラクション0)が成り立つようにする必要がある。 waterPMLコマンドではこの条件を以下の方法で達成する。
The free surface boundary condition (zero traction) is applied on the ground surface (a boundary between solid and vacuum cells) or on the water surface (a boundary between water and vacuum cells). This boundary condition is achieved in the waterPML command by the following way.

4-4-1. 応力の非対角成分 (Off-diagonal stress components)

応力の非対角成分は4つの格子セルが接する辺の中心点で定義されており、地表面や水面が含まれる。 これらの地点において応力が初期値(0)から更新されないようにすることで 自由表面境界条件を実現する(図2)。
Off-diagonal stress components are difined at the center of each edge where four grid cells contact, which includes the ground and water surfaces. The free surface boundary condition is achieved by skipping updates for the off-diagonal stress components from their initial value of zero at these locations (Fig. 2).
プログラム内での具体的な実装方法としては、 応力の定義点のみの通し番号igt(方法3-1) と 半格子点の通し番号ig(方法2) の対応表(igt2ig)を定義する際に 固体内、水中、水底(固体セルと水セルの境界)のみに igtを割り当てるようにしている。
The program implements this approach within the definition of igt2ig, which relates a consecutive index igt for stress definition points (method 3-1) and a consecutive index ig for half-grid nodes (method 2); the values of igt are assigned only to the locations in the solid, in the water, or at the water floor (a boundary between solid and water cells).


図2. 応力の非対角成分を計算する地点 (▪)、 計算をスキップして初期値(0)のままにする地点 (▫)。
Fig. 2. Locations where the off-diagonal stress components are computed (▪) and where the computation is skipped to keep them at their initial values of zero (▫).

4-4-2. 応力の対角成分 (Diagonal stress components)

応力の対角成分は格子セルの中心点で定義されるため、 地表面や水面において応力の対角成分を直接0にすることができない。 そこで、地表面や水面を挟んだ反対側(真空内)の格子セル中心点での法線応力が 固体・水側の格子セル中心点での法線応力の\(-1\)倍になるようにする(図3)。
Because diagonal stress components are defined at the centers of grid cells, they cannot be directly set as zero on the ground and water surfaces. Therefore, the waterPML command assumes that the diagonal stress component at the center of the grid cell that is in the opposite side of the ground or water surface (i.e., the cell in the vacuum) is \(-1\) times that at the center of the solid or water cell (Fig. 3).
実際には真空セル中心点での法線応力を直接設定するのではなく、 地表面や水面における速度の法線方向成分の計算時に間接的にこの条件を与える。 2節で導出した(9)式 \[\begin{eqnarray} V_i^k(\posx_{i-1/2},t) &=& c_0^k(\posx_{i-1/2})V_i^k(\posx_{i-1/2},t_{-1}) \nonumber \\ & & +c_1^k(\posx_{i-1/2})\left[ \sum_l \left\{\tau_{ik}^l\left(\posx_{i-1/2}^{k+1/2},t_{-1/2}\right) -\tau_{ik}^l\left(\posx_{i-1/2}^{k-1/2},t_{-1/2}\right)\right\} \right. \nonumber \\ & & \left. +\delta_{k0}|\Delta\posx_k|f_i(\posx_{i-1/2},t_{-1/2}) \right] \label{eq.Vik} \end{eqnarray}\] において\(\posx_{i-1/2}\)が地表面または水面に位置する場合、 \[\begin{equation} \tau_{ii}^l\left(\posx_{i-1/2}^{i+1/2},t_{-1/2}\right)= -\tau_{ii}^l\left(\posx_{i-1/2}^{i-1/2},t_{-1/2}\right) \label{eq.boundary.normal} \end{equation}\] とおくことで \[\begin{equation} \tau_{ii}^l\left(\posx_{i-1/2}^{i+1/2},t_{-1/2}\right) -\tau_{ii}^l\left(\posx_{i-1/2}^{i-1/2},t_{-1/2}\right) =2\tau_{ii}^l\left(\posx_{i-1/2}^{i+1/2},t_{-1/2}\right) =2\tau_{ii}^l\left(\posx_i,t_{-1/2}\right) \label{eq.boundary.normal.apply.positive_solid} \end{equation}\] \[\begin{equation} \tau_{ii}^l\left(\posx_{i-1/2}^{i+1/2},t_{-1/2}\right) -\tau_{ii}^l\left(\posx_{i-1/2}^{i-1/2},t_{-1/2}\right) =-2\tau_{ii}^l\left(\posx_{i-1/2}^{i-1/2},t_{-1/2}\right) =-2\tau_{ii}^l\left(\posx_{i-1},t_{-1/2}\right) \label{eq.boundary.normal.apply.negative_solid} \end{equation}\] という関係式が得られる。 境界面の\(+\)側(\(\posx_i\))が固体や水のセルである場合、 (\ref{eq.boundary.normal.apply.positive_solid})式を (\ref{eq.Vik})式に代入して得られる \[\begin{eqnarray} V_i^i(\posx_{i-1/2},t) &=& c_0^i(\posx_{i-1/2})V_i^i(\posx_{i-1/2},t_{-1}) \nonumber \\ & & 2c_1^i(\posx_{i-1/2}) \sum_l \tau_{ii}^l\left(\posx_i,t_{-1/2}\right) \nonumber \\ & & \left. +\delta_{i0}|\Delta\posx_i|f_i(\posx_{i-1/2},t_{-1/2}) \right] \label{eq.Vii.boundary.positive_solid} \end{eqnarray}\] という式を用いて速度場を更新する。また境界面の\(-\)側が固体や水のセルである場合、 (\ref{eq.boundary.normal.apply.negative_solid})式を (\ref{eq.Vik})式に代入して得られる \[\begin{eqnarray} V_i^i(\posx_{i-1/2},t) &=& c_0^i(\posx_{i-1/2})V_i^i(\posx_{i-1/2},t_{-1}) \nonumber \\ & & -2c_1^i(\posx_{i-1/2}) \sum_l \tau_{ii}^l\left(\posx_{i-1},t_{-1/2}\right) \nonumber \\ & & \left. +\delta_{i0}|\Delta\posx_i|f_i(\posx_{i-1/2},t_{-1/2}) \right] \label{eq.Vii.boundary.negative_solid} \end{eqnarray}\] という式を用いて速度場を計算する。 いずれの場合にも固体・水側のセルの応力のみを用いて速度場の更新が可能であり、 これらの式は図3の状況を表しているので自由表面境界条件が間接的に満たされる。 このアプローチはOhminato and Chouet (1997)を踏襲したものである。
In practice, the diagonal stress components at the centers of vacuum cells are not directly given; instead, this condition is indirectly used in the computation of the normal velocity component as follows. Suppose that a position \(\posx_{i-1/2}\) is on the ground or water surface in Eq. (9) of section 2 (reproduced as Eq. \ref{eq.Vik}). Assuming that normal stress components around this location satisfy Eq. (\ref{eq.boundary.normal}), we obtain Eqs. (\ref{eq.boundary.normal.apply.positive_solid}) and (\ref{eq.boundary.normal.apply.negative_solid}). If the positive side of the boundary (\(\posx_i\)) is a solid or water cell, the velocity field is updated using Eq. (\ref{eq.Vii.boundary.positive_solid}), which is obtained by inserting Eq. (\ref{eq.boundary.normal.apply.positive_solid}) into (\ref{eq.Vik}). If the negative side of the boundary (\(\posx_{i-1}\)) is a solid or water cell, the velocity field is updated using Eq. (\ref{eq.Vii.boundary.negative_solid}), which is obtained by inserting Eq. (\ref{eq.boundary.normal.apply.negative_solid}) into (\ref{eq.Vik}). In either cases, the velocity field can be updated using the stresses only in the solid or water cells By this, the free-surface boundary condition is satisfied indirectly, as Fig. 3 shows. This approach was introduced by Ohminato and Chouet (1997).


図3. 地表面(や水面)で法線応力を0にするための考え方。 ▪ が計算に登場する固体セル内の法線応力、 ▫ が仮想的に考える真空中の法線応力であり、 後者が前者の\(-1\)倍であるものとして 境界面での速度の法線成分(青)を計算する。
Fig. 3. A schematic illustration for the free surface boundary condition for a normal stress component on a ground (or water) surface. An imaginary normal stress in the vacuum (▫) is assumed to be \(-1\) times that in the solid cell (▪) in the computation of the normal velocity component on the boundary (blue).

このアプローチでは尾根上やピーク上の観測点(図4の赤枠で示す格子セルの中心点)において 速度\(V_i\)が0になってしまうという課題がある。 このことは次のようにして分かる。 この観測点での速度\(V_i\)の計算式は \[\begin{eqnarray} V_i(\posx_i,t) &=& \frac{1}{2}\left[V_i(\posx_{i+1/2},t)+V_i(\posx_{i-1/2},t)\right] \nonumber \\ &=& \frac{1}{2}\sum_k \left[V_i^k(\posx_{i+1/2},t)+V_i^k(\posx_{i-1/2},t)\right] \nonumber \\ &=& \frac{1}{2}\left[V_i^i(\posx_{i+1/2},t)+V_i^i(\posx_{i-1/2},t)\right] +\frac{1}{2}\sum_{k\neq i} \left[V_i^k(\posx_{i+1/2},t)+V_i^k(\posx_{i-1/2},t)\right] \label{eq.Vi.ridge} \end{eqnarray}\] である。 観測点と力源が同じ場所あるいは隣接している特殊なケースを除けば \(f_i(\posx_{i\pm 1/2},t_{-1/2})=0\) であるので(\ref{eq.Vik})式より \[\begin{eqnarray} V_i^k(\posx_{i\pm 1/2},t) &=& c_0^k(\posx_{i\pm 1/2})V_i^k(\posx_{i-1/2},t_{-1}) \nonumber \\ & & +c_1^k(\posx_{i\pm 1/2})\left[ \sum_l \left\{\tau_{ik}^l\left(\posx_{i\pm 1/2}^{k+1/2},t_{-1/2}\right) -\tau_{ik}^l\left(\posx_{i\pm 1/2}^{k-1/2},t_{-1/2}\right)\right\} \right] \label{eq.Vik.no_force} \end{eqnarray}\] が成り立つ。\(i\neq k\)のとき \(\posx_{i\pm 1/2}^{k\pm 1/2}\)は地表面上の点であるので \(\tau_{ik}^l\left(\posx_{i\pm 1/2}^{k\pm 1/2},t_{-1/2}\right)=0\)であり、 したがって速度の初期値が0であれば(\ref{eq.Vik.no_force})式より \(V_i^k(\posx_{i\pm 1/2},t)=0\)となる。 これより(\ref{eq.Vi.ridge})式の右辺第2項は0である。また (\ref{eq.Vii.boundary.positive_solid})(\ref{eq.Vii.boundary.negative_solid}) 式を今回のケースに当てはめると \[\begin{eqnarray} V_i^i(\posx_{i+1/2},t) &=& c_0^i(\posx_{i+1/2})V_i^i(\posx_{i+1/2},t_{-1}) \nonumber \\ & & -2c_1^i(\posx_{i+1/2}) \sum_l \tau_{ii}^l\left(\posx_i,t_{-1/2}\right) \label{eq.Vii.ridge.positive} \end{eqnarray}\] \[\begin{eqnarray} V_i^i(\posx_{i-1/2},t) &=& c_0^i(\posx_{i-1/2})V_i^i(\posx_{i-1/2},t_{-1}) \nonumber \\ & & 2c_1^i(\posx_{i-1/2}) \sum_l \tau_{ii}^l\left(\posx_i,t_{-1/2}\right) \label{eq.Vii.ridge.negative} \end{eqnarray}\] となるので、物性値が同じ (\(c_0^i(\posx_{i+1/2})=c_0^i(\posx_{i-1/2})\), \(c_1^i(\posx_{i+1/2})=c_1^i(\posx_{i-1/2})\)) であれば\(V_i^i(\posx_{i+1/2},t)=-V_i^i(\posx_{i-1/2},t)\)となり、 (\ref{eq.Vi.ridge})式の右辺第1項も0である。 以上より\(V_i(\posx_i,t)=0\)となってしまう。 このプログラム(ymaeda_opentools version 2025-03-25a以降)では このような位置に観測点を置いて波形を計算しようとした場合には 計算をエラー終了する仕様にしている。
This approach results in zero velocity at a station on a ridge or a peak (the center of a grid cell shown by the red frame in Fig. 4). The equation for \(V_i\) at this location is given by Eq. (\ref{eq.Vi.ridge}). The body force \(f_i(\posx_{i\pm 1/2},t_{-1/2})=0\) and thus Eq. (\ref{eq.Vik}) reduces to (\ref{eq.Vik.no_force}) except for limited cases where the source and station are extremely near. The stress \(\tau_{ik}^l\left(\posx_{i\pm 1/2}^{k\pm 1/2},t_{-1/2}\right)=0\) for \(i\neq k\) because \(\posx_{i\pm 1/2}^{k\pm 1/2}\) is on a ground surface. Therefore, Eq. (\ref{eq.Vik.no_force}) indicates that \(V_i^k(\posx_{i\pm 1/2},t)=0\), and thus the 2nd term of the right hand side of Eq. (\ref{eq.Vi.ridge}) is zero. Applying Eqs. (\ref{eq.Vii.boundary.positive_solid}) and (\ref{eq.Vii.boundary.negative_solid}) to this case results in (\ref{eq.Vii.ridge.positive}) and (\ref{eq.Vii.ridge.negative}), and thus \(V_i^i(\posx_{i+1/2},t)=-V_i^i(\posx_{i-1/2},t)\) as long as the material properties at the two locations are same (\(c_0^i(\posx_{i+1/2})=c_0^i(\posx_{i-1/2})\), \(c_1^i(\posx_{i+1/2})=c_1^i(\posx_{i-1/2})\)). Therefore, the 1st term of the right hand side of Eq. (\ref{eq.Vi.ridge}) is also zero. From these, \(V_i(\posx_i,t)=0\) always holds at the stations on a ridge or a peak. This program treats these cases as an error in ymaeda_opentools versions 2025-03-25a and later.


図4. 速度の水平成分が必ず0になってしまう格子セル(赤枠)。
Fig. 4. Grid cells (shown by the red frame) where the horizontal velocity is necessarily zero.