数値モデリング（地質学）

オイラー的アプローチでは、固定された位置における質量や速度などの物理量の時間変化を考慮します。 ^{[ 14 ]}これは、川の水が橋を通過する様子を観察することに似ています。数学的には、物理​​量は位置と時間の関数として表現できます。このアプローチは、自然な境界を持たない流体や均質（均一）な物質に有効です。^{[ 16 ]}

一方、ラグランジュ的アプローチでは、固定された物質要素の体積などの物理量の時間的変化を考慮します。^{[ 14 ]}これは、川を下流に流れる特定の水分子の集合を観察することに似ています。ラグランジュ的アプローチを使用すると、周囲から分離する自然な境界を持つ固体物体を追跡しやすくなります。^{[ 16 ]}

連続方程式は、地質学的対象物または媒体が連続的であることを示す数学的な表現であり、対象物には空隙が存在しないことを意味します。^{[ 17 ]}この方程式は、地質学における数値モデリングでよく使用されます。^{[ 17 ]}

一例として、流体の質量の連続方程式がある。質量保存の法則に基づくと、密度を持つ流体の場合、${\displaystyle \rho }$ポジション${\displaystyle x_{j}}$一定体積内で${\displaystyle V}$流体の場合、質量の変化率は境界を横切る外向きの流体流量に等しい。${\displaystyle S}$:

$\displaystyle {\frac {\partial }{\partial t}}\int \limits _{V}\rho dV=-\int \limits _{S}\rho u_{j}dS_{j}}$

どこ${\displaystyle V}$は体積要素であり、$u_j$速度は${\displaystyle x_{j}}$。

ラグランジュ形式では：^{[ 2 ]}

$${\displaystyle {\frac {D\rho }{Dt}}\equiv {\frac {\partial \rho }{\partial t}}+u_{j}{\frac {\partial \rho }{\partial x_{j}}}=-\rho {\frac {\partial u_{j}}{\partial x_{j}}}}$$

オイラー形式では：^{[ 2 ]}

${\displaystyle {\frac {\partial \rho }{\partial t}}={\frac {\partial (\rho u_{j})}{\partial x_{j}}}}$

この方程式は、モデルが地質学的時間スケールにおけるマントルのような連続的な流体流動を伴う場合に有用である。^{[ 2 ]}

運動量方程式は、加えられた力に応じて物質がどのように動くかを記述する。これはニュートンの運動第二法則の表現である。^{[ 17 ]}

一定の体積を考えてみましょう${\displaystyle V}$物質の。運動量保存の法則により、体積の変化率は次の式に等しい。^{[ 2 ]}

• 外力${\displaystyle F}$要素に適用
• 表面に加わる垂直応力とせん断応力${\displaystyle S}$要素を境界付ける
• その表面上の要素から運動量が減少する。

$${\displaystyle {\frac {\partial }{\partial t}}\int \limits _{V}\rho u_{i}dV=\int \limits _{V}\rho F_{i}dV+\int \limits _{S}\sigma _{ij}dS_{j}-\int \limits _{S}\rho u_{i}u_{j}dS_{j}}$$

どこ${\displaystyle V}$は体積要素です。${\displaystyle u}$は速度です。

簡略化と統合後、あらゆる量に対して${\displaystyle V}$この方程式のオイラー形式は次のとおりです。^{[ 2 ] [ 17 ]}

$\displaystyle \rho {\frac {\partial u_{i}}{\partial t}}+\rho u_{j}{\frac {\partial u_{i}}{\partial x_{j}}}=\rho F_{i}+{\frac {\partial \sigma _{ij}}{\partial x_{j}}}}$

熱方程式は、システム内における熱エネルギーの流れ方を記述する。

エネルギー保存の法則から、エネルギーの変化率は${\displaystyle E}$一定体積の${\displaystyle V}$質量は以下と等しい：^{[ 2 ]}

• 境界で行われた作業${\displaystyle S}$
• 外部の力によって行われた仕事${\displaystyle F}$巻${\displaystyle V}$
• 境界を横切る熱伝導を差し引いた値${\displaystyle S}$
• 境界を横切る熱対流が減少${\displaystyle S}$
• 加えて、体内で発生する熱

数学的に言うと：

$\left\{\displaystyle {\frac {\partial }{\partial t}}\int \limits _{V}\rho EdV=\int \limits _{S}u_{i}\sigma _{ij}dS_{j}+\int \limits _{V}\rho u_{i}F_{i}dV-\int \limits _{S}k{\frac {\partial T}{\partial x_{j}}}dS_{j}-\int \limits _{S}\rho Eu_{j}dS_{j}+\int \limits _{V}\rho HdV}$

どこ${\displaystyle V}$は体積要素です。${\displaystyle u}$速度は、 ${\displaystyle T}$温度は、 ${\displaystyle k}$は伝導係数であり、${\displaystyle H}$は熱発生率である。^{[ 2 ]}

有限要素法では、対象物（または領域）をより小さく、重なり合わない要素（または部分領域）に分割し、これらの要素を節点で接続します。偏微分方程式の解は、より単純な要素方程式（通常は多項式）で近似されます。^{[ 2 ] [ 20 ] [ 21 ]}次に、これらの要素方程式を組み合わせて対象物全体の方程式を作成します。つまり、各要素の寄与を合計して対象物全体の応答をモデル化します。^{[ 2 ] [ 20 ] [ 21 ]}この方法は、機械的な問題を解くためによく使用されます。^{[ 21 ]}有限要素法を使用する一般的な手順は次のとおりです。^{[ 21 ]}

1. 要素タイプを選択してオブジェクトを分割します。一般的な要素タイプには、三角形、四角形、四面体などがあります。^{[ 21 ]}異なる問題には、異なるタイプの要素を選択する必要があります。
2. 変位関数を決定します。変位関数は、要素がどのように動くかを決定します。線形、二次、または三次多項式関数が一般的に使用されます。^{[ 21 ]}
3. 変位とひずみの関係を決定します。要素の変位は、技術的にはひずみと呼ばれる要素の形状の変化または変形を引き起こします。この関係は、変位によって要素がどれだけのひずみを受けたかを計算します。^{[ 21 ]}
4. ひずみ-応力関係を決定します。要素の変形は要素に応力を生じさせます。応力とは、要素に加わる力のことです。この関係は、ひずみによって要素が受ける応力の量を計算します。この関係の例の1つはフックの法則です。^{[ 21 ]}
5. 要素の剛性と剛性マトリックスの式を導出します。応力によって要素が変形します。要素の剛性（剛性）は、応力に応じてどれだけ変形するかを示します。要素のさまざまな方向の剛性は、計算時の操作を簡素化するためにマトリックス形式で表されます。^{[ 21 ]}
6. 要素方程式を組み合わせて全体方程式を作成します。各要素の寄与を合計して、システム全体を記述する一連の方程式を作成します。^{[ 21 ]}
7. 境界条件を適用する。温度、応力、その他の物理量などの境界における事前定義された条件がシステムの境界に導入される。^{[ 21 ]}
8. 変位を解きます。時間が経過するにつれて、要素の変位が段階的に解かれます。^{[ 21 ]}
9. ひずみと応力を求めます。変位が計算された後、ステップ 3 と 4 の関係を使用してひずみと応力を計算します。^{[ 21 ]}

スペクトル法は有限要素法に似ています。^{[ 22 ] [ 23 ]}主な違いは、スペクトル法では基底関数を使用する点です。これは、多数の単純な関数の和によって関数を近似する高速フーリエ変換（FFT）を使用する場合があります。 ^{[ 22 ] [ 23 ]}これらの種類の基底関数は、領域全体に適用され、支配偏微分方程式を近似することができます。^{[ 2 ] [ 22 ] [ 23 ]}したがって、各計算では領域全体の情報が考慮されますが、有限要素法では近傍の情報のみが考慮されます。^{[ 22 ] [ 23 ]}その結果、スペクトル法は指数関数的に収束し、時間または空間の変動が大きい問題の解決に適しています。^{[ 22 ] [ 23 ]}

有限体積法も有限要素法と似ています。対象物をより小さな体積（または要素）に分割し、制御体積上で物理量を、異なる面を通過するこれらの量のフラックスとして解きます。^{[ 2 ] [ 24 ]}使用される方程式は通常、質量やエネルギーなどの物理量の保存またはバランスに基づいています。^{[ 24 ] [ 25 ]}

有限体積法は、有限要素法と同様に不規則なメッシュに適用できます。要素方程式は物理的に意味を持ちます。ただし、要素方程式の高次バージョンは明確に定義されていないため、より高い精度を得ることは困難です。^{[ 2 ] [ 24 ] [ 25 ]}

有限差分法は、微分を差分方程式で近似することにより微分方程式を近似するものであり、偏微分方程式を解く主要な方法である。^{[ 26 ] [ 27 ] [ 28 ] [ 29 ]}

関数を考える${\displaystyle f(x)}$連続かつ有限な関数である一価導関数を持つ${\displaystyle x}$テイラーの定理によれば：^{[ 30 ]}

${\displaystyle f(x+\Delta x)=f(x)+\Delta xf'(x)+{\frac {1}{2}}\Delta x^{2}f''(x)+{\frac {1}{6}}\Delta x^{3}f'''(x)+\cdots }$

そして

${\displaystyle f(x-\Delta x)=f(x)-\Delta xf'(x)+{\frac {1}{2}}\Delta x^{2}f''(x)-{\frac {1}{6}}\Delta x^{3}f'''(x)+\cdots }$

上記の表現をまとめると次のようになります。^{[ 30 ]}

${\displaystyle f(x+\Delta x)+f(x-\Delta x)=2f(x)+\Delta x^{2}f''(x)+{\text{\Delta x の 4 乗を超える項}$

4 乗より大きい項は無視する${\displaystyle x}$すると：^{[ 30 ]}

${\displaystyle f''(x)\simeq {\frac {1}{\Delta x^{2}}}\left[f(x+\Delta x)-2f(x)-f(x-\Delta x)\right]}$

${\displaystyle f'(x)\simeq {\frac {1}{2\Delta x}}\left[f(x+\Delta x)-f(x-\Delta x)\right]}$

上記は導関数の中心差分近似であり、 ^{[ 30 ]}これは前方差分によっても近似できる。

${\displaystyle f'(x)\simeq {\frac {1}{\Delta x}}\left[f(x+\Delta x)-f(x)\right]}$

または後方差分：

${\displaystyle f'(x)\simeq {\frac {1}{\Delta x}}\left[f(x)-f(x-\Delta x)\right]}$

有限差分法の精度は、より高次の項を用いることで向上させることができる。

離散要素法（個別要素法とも呼ばれる）は、不連続性の特性を明示的にモデル化できるため、岩石の節理や層理などの亀裂などの不連続材料をモデル化するためによく使用されます。^{[ 19 ]}この方法は、当初は岩盤力学の問題をシミュレートするために開発されました。 ^{[ 19 ] [ 31 ]}

この方法の主な考え方は、対象物をより小さな粒子の集合体としてモデル化することです。^{[ 19 ]}これは砂で城を建てることに似ています。これらの粒子は球体のような単純な形状をしています。速度などの各粒子の物理量は、粒子間の接触で継続的に更新されます。^{[ 19 ]}このモデルは、多数の粒子を使用する必要があるため、計算負荷が比較的高くなっています。 ^{[ 19 ]}特に斜面のような大規模モデルではその傾向が顕著です。 ^{[ 32 ]}したがって、このモデルは通常、小規模な対象物に適用されます。

数値モデリングは岩盤力学のさまざまな分野で広く応用されている。^{[ 3 ]}岩石は通常次のような性質を持つため、モデル化が難しい材料である。^{[ 3 ]}

• 不連続性：岩盤には多数の亀裂や微細亀裂があり^{[ 37 ]}、岩盤内の空間は空気や水などの他の物質で満たされている可能性があります。^{[ 3 ]}不連続性は岩盤に大きな影響を与えるため、これらの不連続性を完全に捉えるには複雑なモデルが必要です。^{[ 3 ]}
• 異方性：岩盤の特性、例えば透水性（流体が通過できる能力）は、方向によって変化する可能性がある。^{[ 3 ] [ 37 ]}
• 不均質：岩塊の異なる部分の特性が異なる場合がある。^{[ 3 ] [ 37 ]}例えば、花崗岩では石英粒子と長石粒子の物理的特性が異なる。^{[ 38 ] [ 39 ]}
• 非弾性：岩石は応力が取り除かれた後、完全に元の形状に戻ることはできない。^{[ 37 ] [ 3 ]}

岩石の挙動をモデル化するには、上記の特性をすべて考慮に入れた複雑なモデルが必要です。^{[ 3 ]}有限差分法、有限要素法、境界要素法などの方法を使用して岩石を連続体としてモデル化するモデルは多数あります。欠点の1つは、これらのモデルでは通常、亀裂やその他の不連続性をモデル化する能力が限られていることです。^{[ 40 ]}離散要素法や離散破壊ネットワーク法などの方法を使用して岩石を不連続体としてモデル化するモデルも一般的に使用されています。^{[ 3 ] [ 35 ]}両方の方法の組み合わせも開発されています。^{[ 3 ]}

数値モデリングは、数値実験を行うことで岩石の力学的プロセスの理解を深め、設計や建設作業に役立ちます。^{[ 3 ]}

数値モデリングは、地球の地殻の熱履歴を予測および記述するために使用されており、地質学者が熱年代学的データの解釈を改善できる。^{[ 41 ]}熱年代学は、岩石が特定の温度以下に冷却された時間を示すことができる。^{[ 42 ]}断層の形成や表面侵食などの地質学的イベントは、表面で採取されたサンプルの熱年代学的パターンを変化させる可能性があり、これらのデータによって地質学的イベントを制約することができる。^{[ 42 ]}数値モデリングは、パターンを予測するために使用できる。

地殻の熱モデリングの難しさは、主に地表の不規則性と時間経過に伴う変化（主に侵食）に関係しています。したがって、地表の形態変化をモデル化するために、モデルは時間とともに変化する境界条件と不規則なメッシュを持つ熱方程式を解く必要があります。 ^{[ 43 ]}

ペキュベ

Pecubeは、熱年代学的パターンを予測するために開発された数値モデルの1つです。^{[ 43 ]}これは、有限要素法を使用して、移流を伴う以下の一般化された熱伝達方程式を解きます。^{[ 41 ]}右辺の最初の3つの項は、伝導によって伝達される熱です。${\displaystyle x}$、${\displaystyle y}$そして${\displaystyle z}$指示に従って${\displaystyle A}$移流です。 $${\displaystyle {\frac {\partial T}{\partial t}}+u{\frac {\partial T}{\partial x}}+v{\frac {\partial T}{\partial y}}+w{\frac {\partial T}{\partial z}}={\frac {\partial }{\partial x}}\kappa {\frac {\partial T}{\partial x}}+{\frac {\partial }{\partial y}}\kappa {\frac {\partial T}{\partial y}}+{\frac {\partial }{\partial z}}\kappa {\frac {\partial T}{\partial z}}+A}$$ モデル内で温度場が構築された後、粒子の経路が追跡され、粒子の冷却履歴が得られます。^{[ 41 ]}その後、熱年代学的年代のパターンを計算できます。^{[ 41 ]}

水文地質学 では、地下水流はしばしば有限要素法^{[ 46 ] [ 47 ] [ 48 ]}と有限差分法^{[ 49 ]}によって数値的にモデル化される。これらの2つの方法は、メッシュが十分に細かい場合、同様の結果を生み出すことが示されている^{[ 50 ] [ 51 ]。}

地殻と岩石圏のレオロジー（応力に対する材料の応答）は、自由表面（地表面）と地殻材料の塑性および弾性を考慮する必要があるため複雑です。[ 2 ]ほとんどのモデルは^{、ラグランジュメッシュ}を使用した有限要素法を使用しています。^{[ 2 ]} 1つの用途は、沈み込みの変形と運動学の研究です。^{[ 54 ] [ 55 ]}

FLAC

連続体の高速ラグランジュ解析 (FLAC)は、地殻ダイナミクスのモデリングで最も一般的なアプローチの 1 つです。^{[ 2 ]}このアプローチは、行列を使用せずに運動量方程式と連続方程式を解くため高速ですが、時間ステップは十分に小さくする必要があります。 ^{[ 56 ]このアプローチは、地殻ダイナミクスの 2D、 [ 57 ] [ 58 ] [ 59 ]} 2.5D、^{[ 60 ]}および 3D ^{[ 61 ]}の研究で使用されており、2.5D の結果は、2 次元の結果の複数のスライスを組み合わせることによって生成されました。^{[ 2 ]}

マントル対流をモデル化しようとする試みは数多く存在する。

有限要素法^{[ 65 ]} 、有限体積法、有限差分法^{[ 66 ]}、スペクトル法はすべてマントル対流のモデリングに使用されており、ほとんどすべてのモデルでオイラー格子が使用されています。^{[ 2 ]}有限差分法とスペクトル法は単純で高速であるため、初期のモデルの一部で使用されていましたが、2010年代には一般的に有限要素法または有限体積法が採用されました。^{[ 2 ]}多くのベンチマーク論文がこれらの数値モデルの妥当性を調査しています。^{[ 2 ] [ 67 ] [ 68 ] [ 69 ] [ 70 ] [ 71 ] [ 72 ]}現在のアプローチでは、ほとんどが固定された均一な格子を使用しています。^{[ 2 ]}より正確な近似が必要な部分で要素のサイズを小さくする格子細分化は、マントル対流の数値モデリングの将来の発展の方向性である可能性があります。^{[ 2 ] [ 73 ]}

スペクトル法

マントル対流におけるスペクトル法は、3 次元の支配方程式をいくつかの 1 次元方程式に分解し、方程式をはるかに速く解く。これは、初期のマントル対流モデルでよく用いられた手法の 1 つです。^{[ 2 ]} 1980 年代から 2000 年代初頭にかけて、この手法を用いた多くのプログラムが開発されました。^{[ 2 ] [ 90 ] [ 91 ] [ 92 ] [ 93 ] [ 94 ] [ 95 ] [ 96 ]}しかし、この手法ではマントルの粘性の横方向の変化を扱うのが難しく、2010 年代には他の手法がより普及しました。^{[ 2 ]}

プレートテクトニクスは、地球の岩石圏が基本的にマントル上に浮かぶプレートで構成されているという理論です。^{[ 97 ]}マントル対流モデルは、その上に浮かぶプレートをモデル化する上で基本であり、プレートをこのモデルに組み込むには、剛体ブロックアプローチとレオロジーアプローチという 2 つの主要なアプローチがあります。^{[ 2 ]}剛体ブロックアプローチでは、プレートは剛体であると仮定します。つまり、プレートは形状を維持し、変形しません。これは、水に浮かぶ木片のようなものです。対照的に、レオロジーアプローチでは、プレートを高粘性流体としてモデル化し、その下の岩石圏に適用される方程式が上のプレートにも適用されます。^{[ 2 ]}

地球の核内の導電性の鉄とニッケルの流体の動きによって地磁気場が生成されるという理論である地磁気ダイナモ理論を検証するために数値モデルが作成されました。^{[ 2 ] [ 98 ]}

地球の液体の外核の流れをモデル化することは困難です。なぜなら、次の理由があるからです。^{[ 2 ]}

• 地球の自転によるコリオリ効果は無視できない
• 発生した磁場はローレンツ力も発生させ、それが液体外核内の導電性流体の運動に影響を与える。
• 液体鉄の粘度が低いため、流体の流れをモデル化するのが困難である。

ほとんどのモデルは、例えばグラッツマイヤー・ロバーツモデルのように、スペクトル法を用いて^ジオダイナモをシミュレートしています。^{[ 2][99} ]また^{、カゲヤマとサトウのモデルで}は、有限差分法も使用されています。^{[ 99 ] [ 102 ]さらに、有限体積法[ 103 ]や有限要素法[ 104 ]など、他の方法を}試した研究もあります。

有限差分法は、地震波の伝播のシミュレーションに広く用いられてきた。^{[ 106 ] [ 107 ] [ 108 ]}しかし、計算能力の制限により、一部のモデルではメッシュの間隔が大きすぎるため（地震波の波長と比較して）、異なる周波数の地震波が分離するグリッド分散により結果が不正確になる。 ^{[ 106 ] [ 109 ]}一部の研究者は、スペクトル法を用いて地震波の伝播をモデル化することを提案している。^{[ 106 ] [ 110 ]}