17:15 〜 18:45
[SGD01-P04] Lagrange形式に基づく地上重力変化の数値計算：球状圧力源変動を例に
キーワード：重力変化、有限要素法、球状圧力源、地殻変動、Lagrange形式
地殻変動観測や地上重力観測といった測地学的観測は、地震火山活動に関連する地球内部変動を把握するのに有効である。これらの観測データを用いて地球内部変動現象をモデル化する際には、断層運動や球状圧力源変動に関する解析解（山川, 1955; 萩原, 1977; Okada, 1992; Okubo, 1992）がグリーン関数としてしばしば用いられる。しかし、これらの解析解は均質等方な半無限媒質を仮定しているため、媒質不均質・地表形状・地球曲率などが無視できない状況でこれらの解析解をグリーン関数に用いると、モデル結果にバイアスが乗る可能性がある。また、萩原（1977）の解析解は球状圧力源変動に伴う重力変化を計算する際に地表面の形状変化を無限平板で近似するため、現実の地上重力変化を十分に再現できていない可能性がある。
上述の問題点を解決しうる方法として有限要素法がある。有限要素法では、地形や地下構造などの複雑性を力学的に取り入れて地殻変動や重力変化を計算することができる。近年では地殻変動をモデル化する際に有限要素法が頻繁に利用されているほか（Iinuma et al., 2016）、地殻変動に関する解析解の適用限界を有限要素法の立場から考察している例もある（坂井他, 2007）。一方、重力変化は、有限要素解析で密度変化場と観測点の鉛直変位を求めた上で、境界値問題の求解（Currenti, 2007）や空間積分（Trasatti and Bonafede, 2008）によって計算されてきた。しかし、これらの手法はEuler形式に基づいて観測点の移動を単純化しているため、実際の地上重力変化を精度よく再現できない可能性がある。また、これらの手法はフリーエア勾配を定数として扱っており、地形や媒質不均質がフリーエア勾配値に与える影響を考慮していない。
そこで本研究は、重力変化計算に関する複数の問題を解決するために、有限要素法を活用して地上重力変化を計算する新たな手法を提案する。本手法では、まず各要素に対して質量保存条件を維持しながら変形前後の密度変化量をLagrange形式で記述する。次に、変形前後の密度場における万有引力効果（g1およびg2）を空間積分によって計算し、それらの差分を取ることで有限要素モデル由来の重力変化(dgFEM = g2 – g1)を計算する。この際、重力観測点は媒質表面の任意の要素上に固定するので、観測点の移動を厳密に扱いつつ、有限要素モデルから受けるフリーエア効果を含めて重力変化を計算することができる。さらに、有限要素解析で得られた観測点の上下変位dhFEMを用いて、有限要素モデルの領域外から受けるフリーエア効果をdgres = dhFEM * βresと計算する。なお、βresはモデル領域外の質量分布がもたらすフリーエア勾配であり、質量分布の幾何学的対称性から別途計算する。以上の計算により、地上重力変化の観測値と比較できる値として、重力変化の総和dg = dgFEM + dgresを得る。
我々は本手法の有効性を評価するために、半無限媒質の条件下における球状圧力源変動の有限要素モデルを用いて地上重力変化を計算し、これを解析解と比較した。なお、地表変位量および空中固定点での重力変化量については、有限要素モデルの数値計算結果が解析解（山川, 1955; 萩原, 1977）と誤差数％以内で一致することを確認済みである。その結果、本手法で計算された地上重力変化についても、有限要素解析の誤差範囲内で萩原（1977）の解析解と一致することが分かった。このことは、萩原（1977）で実施していた地表変形の無限平板近似が、地上重力変化を計算する上でも十分な精度で成り立つことを意味している。ただし、これらは半無限媒質中の球状圧力源変動を扱った場合の結果であり、地表の凹凸を加味した場合には有限要素モデルにおける地上重力変化が解析解と乖離する可能性が高い。そこで本研究は、今後有限要素モデルに円錐地形を採用した上で地上重力変化を計算し、それを解析解（Nishiyama, 2022）と比較する予定である。
上述の問題点を解決しうる方法として有限要素法がある。有限要素法では、地形や地下構造などの複雑性を力学的に取り入れて地殻変動や重力変化を計算することができる。近年では地殻変動をモデル化する際に有限要素法が頻繁に利用されているほか（Iinuma et al., 2016）、地殻変動に関する解析解の適用限界を有限要素法の立場から考察している例もある（坂井他, 2007）。一方、重力変化は、有限要素解析で密度変化場と観測点の鉛直変位を求めた上で、境界値問題の求解（Currenti, 2007）や空間積分（Trasatti and Bonafede, 2008）によって計算されてきた。しかし、これらの手法はEuler形式に基づいて観測点の移動を単純化しているため、実際の地上重力変化を精度よく再現できない可能性がある。また、これらの手法はフリーエア勾配を定数として扱っており、地形や媒質不均質がフリーエア勾配値に与える影響を考慮していない。
そこで本研究は、重力変化計算に関する複数の問題を解決するために、有限要素法を活用して地上重力変化を計算する新たな手法を提案する。本手法では、まず各要素に対して質量保存条件を維持しながら変形前後の密度変化量をLagrange形式で記述する。次に、変形前後の密度場における万有引力効果（g1およびg2）を空間積分によって計算し、それらの差分を取ることで有限要素モデル由来の重力変化(dgFEM = g2 – g1)を計算する。この際、重力観測点は媒質表面の任意の要素上に固定するので、観測点の移動を厳密に扱いつつ、有限要素モデルから受けるフリーエア効果を含めて重力変化を計算することができる。さらに、有限要素解析で得られた観測点の上下変位dhFEMを用いて、有限要素モデルの領域外から受けるフリーエア効果をdgres = dhFEM * βresと計算する。なお、βresはモデル領域外の質量分布がもたらすフリーエア勾配であり、質量分布の幾何学的対称性から別途計算する。以上の計算により、地上重力変化の観測値と比較できる値として、重力変化の総和dg = dgFEM + dgresを得る。
我々は本手法の有効性を評価するために、半無限媒質の条件下における球状圧力源変動の有限要素モデルを用いて地上重力変化を計算し、これを解析解と比較した。なお、地表変位量および空中固定点での重力変化量については、有限要素モデルの数値計算結果が解析解（山川, 1955; 萩原, 1977）と誤差数％以内で一致することを確認済みである。その結果、本手法で計算された地上重力変化についても、有限要素解析の誤差範囲内で萩原（1977）の解析解と一致することが分かった。このことは、萩原（1977）で実施していた地表変形の無限平板近似が、地上重力変化を計算する上でも十分な精度で成り立つことを意味している。ただし、これらは半無限媒質中の球状圧力源変動を扱った場合の結果であり、地表の凹凸を加味した場合には有限要素モデルにおける地上重力変化が解析解と乖離する可能性が高い。そこで本研究は、今後有限要素モデルに円錐地形を採用した上で地上重力変化を計算し、それを解析解（Nishiyama, 2022）と比較する予定である。