地球動力學(xué)數(shù)值模擬算法的應(yīng)用現(xiàn)狀與展望

摘要： 當(dāng)前，數(shù)值模擬技術(shù)已成為探究地球動力學(xué)過程的重要手段。本文綜述了數(shù)值模擬在地球動力學(xué)研究中的應(yīng)用現(xiàn)狀、主要算法和代碼并探討了發(fā)展前景。數(shù)值模擬通過對連續(xù)方程組進(jìn)行離散化，并利用數(shù)值算法求解，從而高效模擬地質(zhì)過程，預(yù)測地球系統(tǒng)的動態(tài)變化。本文系統(tǒng)闡述了有限元法、有限差分法、邊界元法和有限體積法等主流數(shù)值方法的原理、特點和適用條件；并重點探討了有限元法在構(gòu)造應(yīng)力場和斷層動力學(xué)模擬中的應(yīng)用、有限差分法在地震波傳播模擬中的應(yīng)用、邊界元法和有限體積法在斷層力學(xué)和孔隙流體流動模擬中的應(yīng)用。通過對比分析不同方法的優(yōu)缺點，揭示了多種數(shù)值方法耦合是未來地球動力學(xué)數(shù)值模擬的重要發(fā)展方向；隨著高性能計算技術(shù)的進(jìn)步和大數(shù)據(jù)時代的到來，數(shù)值模擬技術(shù)將在地球系統(tǒng)科學(xué)研究中發(fā)揮越來越重要的作用。

關(guān)鍵詞：數(shù)值模擬；地球動力學(xué)；地質(zhì)建模；數(shù)值方法；地球系統(tǒng)；有限元法；多物理場耦合模擬

doi：10.13278/j.cnki.jjuese.20240106 中圖分類號：P541 文獻(xiàn)標(biāo)志碼：A

收稿日期： 2024-05-12

作者簡介： 魏虹羽（1996—），女， 碩士研究生， 主要從事礦產(chǎn)普查與勘探方面的研究，E-mail：weihy21@mails.jlu.edu.cn

通信作者： 李世超（1980—）， 男， 教授，博士，博士生導(dǎo)師，主要從事構(gòu)造地質(zhì)學(xué)方面的研究， E-mail： lsc@jlu.edu.cn

基金項目： 國家自然科學(xué)基金項目（41872234）

Supported by the National Natural Science Foundation of China （41872234）

Current Status and Future Prospects of Numerical Simulation Algorithms in Geodynamics

Wei Hongyu ，Li Shichao，Wang Weian

College of Earth Science， Jilin University， Changchun 130061， China

Abstract： Numerical simulation technology has emerged as a crucial tool for exploring geodynamic processes. This paper reviews the current applicationsof numerical simulations in geodynamics， focusing on key methods， algorithms， and codes， with outlining future development prospects. By discretizing continuous equations and employing numerical algorithms， these simulations efficiently model geological processes and predict dynamic changes in the Earth’s system. The paper systematically discusses the principles， characteristics， and applicable conditions for mainstream numerical methods， including the finite element method， finite difference method， boundary element method， and finite volume method. It explores the finite element method for simulating tectonic stress fields and fault dynamics， the finite difference method for seismic wave propagation， and the boundary element and finite volume methods for fault mechanics and pore fluid flow. A comparative analysis highlights the strengths and weaknesses of each method， indicating that coupling various numerical methods is a key future direction in geodynamics. Furthermore， this paper anticipates the growing significance of numerical simulation technology in Earth system scientific research， propelled by advancements in high-performance computing and the era of big data.

Key words： numerical simulation；geodynamic；geological modeling；numerical methods；Earth system；finite element method; multiphysics coupling simulation

0 引言

地球動力學(xué)是研究地球內(nèi)部和表面動力過程及其相互作用的科學(xué)，涉及地質(zhì)學(xué)、地球物理學(xué)、地球化學(xué)等多個學(xué)科領(lǐng)域。近年來，隨著計算機(jī)技術(shù)和數(shù)值模擬方法的快速發(fā)展，地球動力學(xué)數(shù)值模擬已經(jīng)成為認(rèn)識和理解地球動力過程的重要途徑。利用數(shù)值模擬方法可以在不同時空尺度上模擬地球動力過程，為探索地球內(nèi)部結(jié)構(gòu)、物質(zhì)運(yùn)動和能量轉(zhuǎn)換提供了新的視角和途徑。

數(shù)值模擬技術(shù)在地球動力學(xué)多個領(lǐng)域得到了廣泛應(yīng)用。在地質(zhì)學(xué)領(lǐng)域，數(shù)值模擬被應(yīng)用于研究巖石變形、斷層滑移和流體運(yùn)移等問題[1-3]；在地球化學(xué)領(lǐng)域，數(shù)值模擬被應(yīng)用于研究元素遷移、化學(xué)反應(yīng)和流體-巖石相互作用等問題[4-6]；在地球物理學(xué)領(lǐng)域，數(shù)值模擬被應(yīng)用于研究地幔對流、板塊運(yùn)動、地磁場演化、地震波傳播、震源機(jī)制和地震災(zāi)害預(yù)測等問題[7-9]。

盡管地球動力學(xué)數(shù)值模擬如今已經(jīng)取得了巨大進(jìn)展，但仍然面臨著諸多挑戰(zhàn)。首先，地球動力過程涉及從微觀的原子尺度到宏觀的行星尺度等多個時空尺度，如何在不同尺度之間進(jìn)行耦合和跨尺度模擬仍然是一個難題[10]；其次，地球物質(zhì)的材料特性復(fù)雜多變，如何準(zhǔn)確描述巖石圈本構(gòu)關(guān)系和流變學(xué)行為也是一個挑戰(zhàn)[11-12]；此外，如何有效地整合不同來源的觀測數(shù)據(jù)，并用于約束和驗證數(shù)值模型，也是一個亟待解決的問題。但影響數(shù)值模擬方法的有效性和模擬效率的最大挑戰(zhàn)是數(shù)值模擬方法的選擇。

近年來，地球動力學(xué)數(shù)值模擬算法取得了長足進(jìn)展。高性能計算技術(shù)的發(fā)展使得大規(guī)模并行計算成為可能，也極大地提高了數(shù)值模擬的效率和精度[13]。與此同時，各種新的數(shù)值方法也不斷涌現(xiàn)，如自適應(yīng)有限元法、多重邊界元法、Crank-Nicolson方法等，為解決復(fù)雜的地球動力學(xué)問題提供了更多選擇。

本文綜述地球動力學(xué)數(shù)值模擬算法的最新進(jìn)展和應(yīng)用現(xiàn)狀，重點針對地球動力學(xué)數(shù)值模擬中主流數(shù)值計算方法的最新進(jìn)展，闡述其原理、特點和適用條件，以及這些方法在地質(zhì)學(xué)、地球物理學(xué)、地球化學(xué)等領(lǐng)域的應(yīng)用。同時，也簡要介紹了目前流行的成熟模擬軟件/代碼，為地球動力學(xué)數(shù)值模擬的應(yīng)用提供參考。

1 數(shù)值模擬基本概念

數(shù)學(xué)建模通過一組由偏微分方程和特定邊界條件以及初始條件搭建起來的數(shù)學(xué)模型，定量模擬地球動力學(xué)中的復(fù)雜地質(zhì)過程，比如模擬地殼和地幔在地質(zhì)時間尺度上的緩慢變形。這些變形過程通常伴隨著熱傳輸、地球深部的相變、復(fù)雜的流變學(xué)行為、熔融和熔體遷移、化學(xué)反應(yīng)以及固體運(yùn)動等現(xiàn)象。

要求解數(shù)學(xué)模型，首先需要通過前期實驗獲取相應(yīng)的數(shù)據(jù)，設(shè)定模型的初始條件和邊界條件參數(shù)；然后根據(jù)具體問題選擇合適的數(shù)值方法求解相應(yīng)的地球物理過程。雖然解析解可以提供精確的結(jié)果，但由于模型的復(fù)雜性，許多情況下無法通過解析解得到答案。因此，我們使用計算方法和計算機(jī)來獲得近似解。隨著計算機(jī)硬件技術(shù)的不斷進(jìn)步，計算機(jī)底層運(yùn)算方式的改變也促進(jìn)了數(shù)值模擬算法的改進(jìn)，更快、更精確的算法不斷產(chǎn)生，這使得數(shù)值模擬可以解決更為復(fù)雜的地質(zhì)問題。[HJ2mm]

數(shù)值模擬算法使我們能夠利用先進(jìn)的計算機(jī)技術(shù)來模擬和分析現(xiàn)實世界中各種復(fù)雜的地質(zhì)作用過程。搭建數(shù)值模型的步驟如下：①確定物理參數(shù)。收集或測量模擬系統(tǒng)中的關(guān)鍵物理量。②建立物理模型。根據(jù)物理現(xiàn)象構(gòu)建對應(yīng)的物理模型。③建立數(shù)學(xué)模型。將物理模型轉(zhuǎn)化為數(shù)學(xué)模型，通過微分方程描述物理量隨時間和空間的變化，并定義邊界條件和初始條件。④校準(zhǔn)模型準(zhǔn)確性。通過與觀測數(shù)據(jù)的比較，對模型進(jìn)行校正以提高其準(zhǔn)確性。⑤選擇合適的數(shù)值模擬算法。根據(jù)問題的性質(zhì)和需求，選擇合適的數(shù)值方法。⑥構(gòu)建對應(yīng)的數(shù)值代碼。將連續(xù)的數(shù)學(xué)方程轉(zhuǎn)換為離散方程，便于數(shù)值求解。⑦求解。利用計算機(jī)程序?qū)﹄x散方程進(jìn)行求解。⑧可視化。對數(shù)值模擬的結(jié)果進(jìn)行可視化，以便分析和解釋。數(shù)值模擬算法本質(zhì)上是一系列指令。

相對于上述步驟中數(shù)值模擬算法的選擇，我們更關(guān)心模擬結(jié)果的準(zhǔn)確性，需要注意的是，模型的簡化程度、理想化數(shù)學(xué)模型的選擇、求解離散方程時離散質(zhì)量的差異、迭代方法是否充分等因素都有可能導(dǎo)致計算結(jié)果精確度發(fā)生變化。因此，對數(shù)值模擬算法的選擇分析和解釋都需要仔細(xì)考量，比如通過采取更精確的插值近似或?qū)⒔茟?yīng)用于更小的域來減少離散化誤差等。在數(shù)值計算方法的選擇上，有限元法、有限差分法和有限體積法是目前廣泛應(yīng)用于建模領(lǐng)域的主流方法。譜方法在復(fù)雜幾何處理上存在局限性，但也被用于特定情況下的三維建模，比如需要高精度解的中小規(guī)模問題。邊界元法和離散元法則常用于地球動力學(xué)建模，尤其是在處理非均質(zhì)介質(zhì)和復(fù)雜邊界條件時有較大優(yōu)勢。從理論層面考慮，當(dāng)計算網(wǎng)格足夠精細(xì)時，不同的數(shù)值方法都應(yīng)當(dāng)收斂于相同的解。但在實際運(yùn)用過程中，由于計算資源和時間的限制，研究者仍需要根據(jù)具體的研究問題和條件選擇最為合適的數(shù)值方法。這就要求研究者對各種數(shù)值方法的理論基礎(chǔ)、適用條件以及計算效率有深入地理解。

2 理論基礎(chǔ)

2.1 控制方程

地球動力學(xué)中建立模型一般基于質(zhì)量、動量和熱量三大守恒方程，根據(jù)這三大守恒方程揭示速度、壓力和溫度等變量的值，以及它們?nèi)绾坞S時間演變和如何在空間中變化。由于地球動力學(xué)領(lǐng)域研究問題的時間尺度通常在數(shù)千年到數(shù)百萬年，所以建模過程中通常將地球內(nèi)部視為黏性流體，采用斯托克斯方程描述高黏性流體在引力場驅(qū)動下的運(yùn)動。

質(zhì)量守恒方程：

式中：ρ為密度；t為時間；v為速度向量。式（1）表明了任何給定的材料體積內(nèi)流入/流出質(zhì)量隨時間的局部變化。

當(dāng)然根據(jù)研究問題的具體情況，上述質(zhì)量守恒方程也會進(jìn)行相應(yīng)的調(diào)整。如假設(shè)物質(zhì)不可壓縮，或者考慮地球為黏彈性（局部質(zhì)量變化與地球整體質(zhì)量通量對比可以忽略不計）時，方程應(yīng)進(jìn)行簡化來匹配具體的研究問題。

動量守恒方程：

式中：σ為應(yīng)力張量；g為重力加速度矢量。對于地幔對流和長期構(gòu)造模型來說，其他物體力相較于重力可以忽略不計，所以只考慮重力即可。在巖石圈尺度的研究中，g=9.8 m s-2。

式中：Cp為等壓熱容（p為壓力） ；T為溫度；k為導(dǎo)熱系數(shù)；H為體積產(chǎn)熱項；S為其他熱源。熱能量隨時間的變化可以由幾個不同的過程引起。ρ CpT/t項對應(yīng)的是對流熱輸運(yùn)，即熱能隨物質(zhì)運(yùn)動的輸運(yùn)。它取決于物質(zhì)運(yùn)動的速度v，以及溫度和熱能在空間中的變化（上述三大守恒方程需要針對具體模型的物理過程進(jìn)行選擇，以便得出更符合實際的模擬結(jié)果[14]。

2.2 流變學(xué)

除上述守恒方程外，我們還需要對固體地球?qū)ν饬Φ捻憫?yīng)做出假設(shè)，通過本構(gòu)方程將材料性質(zhì)與守恒方程的變量（溫度等）聯(lián)系起來。其中流變性對巖石圈和地殼的模型至關(guān)重要，它體現(xiàn)了外部施加板塊驅(qū)動力區(qū)域模型中應(yīng)力與變形之間的關(guān)系。不同研究問題需要不同的流變學(xué)性質(zhì)。如巖石的非彈性變形在長時間尺度上的變化過程在模擬中一般視為高黏性流體的非彈性變形，通過Stokes方程來描述。物理模型通?？梢猿浞纸忉尩厍騼?nèi)部與地幔對流和觀測有關(guān)的許多過程。但也有一些與流體無關(guān)的過程，如板塊構(gòu)造、涉及應(yīng)變局部化和強(qiáng)滯后等。這些流變學(xué)超出了黏性流體的基本假設(shè)，包括塑性屈服、應(yīng)變、應(yīng)變速率減弱和硬化。因此如何解決黏性流體的假設(shè)也是建模時需要考慮的問題。常見流變學(xué)過程如下。

地殼和地幔內(nèi)較短時間尺度過程的建模應(yīng)考慮彈性性質(zhì)，低溫條件下巖石會產(chǎn)生脆性或塑性變形，由此導(dǎo)致的剪切帶等研究問題需要考慮黏-彈-塑性特性。其他復(fù)雜地質(zhì)變化過程需要在基礎(chǔ)方程上增加相應(yīng)依賴項，使方程能夠準(zhǔn)確描述物理過程。具體推導(dǎo)過程及解析可參考經(jīng)典地球動力學(xué)教科書[15]。

3 數(shù)值計算方法

在現(xiàn)實生活中，宏觀運(yùn)動和微觀運(yùn)動均基于物質(zhì)和能量的基本性質(zhì)。因此數(shù)值模擬的第一步就是搭建相應(yīng)的物理模型（一組方程）來描述對應(yīng)物理或者化學(xué)變化。對于簡單的物理模型，通常使用解析解和半解析解進(jìn)行求解。然而，在地球動力學(xué)研究領(lǐng)域，所研究的問題往往涉及較為復(fù)雜的時空關(guān)系和多物理場耦合。為了與現(xiàn)實相符合，搭建模型方程時需要根據(jù)需求增加多個依賴關(guān)系式，由此產(chǎn)生的復(fù)雜方程需要選擇合適的數(shù)值模擬算法來求解。

3.1 有限元法

有限元法是用于解決由偏微分方程描述或可表述為泛函最小化問題的方法。研究域通常被表示為有限元素的集合。有限元中近似函數(shù)是根據(jù)所尋找物理場的節(jié)點值來確定的。將一個連續(xù)的物理問題轉(zhuǎn)化為具有未知節(jié)點值的離散有限元問題。

有限元法具有以下特性：

1）高精度近似。有限元法可以通過簡單的近似函數(shù)提供較高的近似值，通過變分原理或加權(quán)余值技術(shù)構(gòu)建的有限元方程，確保了數(shù)學(xué)模型與物理現(xiàn)象的一致性，這種方法的精確度可以通過增加元素數(shù)量或提高近似函數(shù)的階數(shù)來提高。誤差為O（hm+1-s），其中h為元素尺寸，m為近似函數(shù)的多項式階數(shù)，s為平滑指數(shù)。

2）收斂性高。有限元法可系統(tǒng)化為矩陣形式運(yùn)算，適合計算機(jī)編程。當(dāng)前計算機(jī)水平已經(jīng)可以處理多物理場耦合的復(fù)雜情況，隨著有限元法的不斷改良，通過網(wǎng)格加密和細(xì)化策略，已開發(fā)出不需要特殊處理振蕩項和內(nèi)節(jié)點性質(zhì)的針對2D彈性問題的自適應(yīng)有限元法[16]。在特定條件下，有限元解或其導(dǎo)數(shù)在某些特殊點上可能表現(xiàn)出超過標(biāo)準(zhǔn)誤差估計的精確性。這種現(xiàn)象被稱為超收斂性，它在實際計算中已被廣泛觀察到，并從理論上得到了研究和證明[17-18]。

3）靈活性和適應(yīng)性。有限元法在幾何形狀、邊界條件等不同的情況下有很高的靈活性和通用性，可以選用任意形狀的單元，且可以針對不同的物理問題選用不同的單元拓?fù)浣Y(jié)構(gòu)。最簡單單元是兩個節(jié)點的一維單元，同時也可以設(shè)置成三角單元或者四面體。此外，有限元法不限制場函數(shù)滿足的方程類型，延展性強(qiáng)，不會局限于某一特定問題。

有限元法是目前最常見的數(shù)值方法，廣泛應(yīng)用于應(yīng)力分析、熱分析、流體動力學(xué)、電磁學(xué)和多物理場耦合等問題中。

最早17世紀(jì)，Newton等[19]確立了積分運(yùn)算的全局可加性質(zhì)。18世紀(jì)，Gauss引入了加權(quán)殘差法和線性代數(shù)系統(tǒng)的解決方案。Galerkin于1915年提出了一種選擇位移函數(shù)中形函數(shù)的方法，即著名的Galerkin方法，該方法后來被廣泛應(yīng)用于有限元法[20]。Hrennikoff [21]首次將有限元法應(yīng)用于求解彈性力學(xué)問題，探索一維桿中應(yīng)力，這是有限元法應(yīng)用的開端。Courant [22]于1943年首先在定義域內(nèi)局部采用位移函數(shù)來表示未知函數(shù)的方法，盡管當(dāng)時計算機(jī)尚未出現(xiàn)，沒有得到廣泛關(guān)注，但這一方法實質(zhì)上奠定了有限元法的基礎(chǔ)。Turner等[23]于1956年首次提出在二維單元中利用單元剛度矩陣分析構(gòu)建整體結(jié)構(gòu)剛度矩陣的程序。Clough [24]在1960發(fā)表的論文中首次引進(jìn)了“有限元”的概念，并在這篇論文中將有限元法應(yīng)用到了土木工程領(lǐng)域。到20世紀(jì)60年代，有限元法幾乎可被應(yīng)用于所有連續(xù)介質(zhì)力學(xué)所涉及的領(lǐng)域。我國馮康教授[25]《基于變分原理的差分格式》一書的發(fā)表代表著有限元的概念首次在中國提出，奠定了中國有限元法的基礎(chǔ)。隨著越來越多的數(shù)學(xué)家參與到有限元法的研究工作中，有限元法逐步脫離了針對無法解決復(fù)雜工程問題的局限性，取而代之的是統(tǒng)一的觀點和嚴(yán)密的數(shù)學(xué)描述，使這種方法被更多科學(xué)家接受并廣泛應(yīng)用于各種實際問題中。

20世紀(jì)60年代末，有限元法被地學(xué)界采用進(jìn)行構(gòu)造解析，并因為其獨特優(yōu)勢被廣泛應(yīng)用于地學(xué)的各項領(lǐng)域中。70年代中期開始，中國地學(xué)界也開始采用有限元數(shù)值模擬技術(shù)，使得科學(xué)家能夠創(chuàng)建精細(xì)的地質(zhì)構(gòu)造模型，以便更直觀地洞察構(gòu)造應(yīng)力場的動態(tài)特性及演化路徑，并對未來的演變趨勢進(jìn)行科學(xué)預(yù)測。尤其隨著計算機(jī)處理速度顯著地提升，數(shù)值模擬技術(shù)實現(xiàn)了跨越式發(fā)展。各項領(lǐng)域的科研工作者可以執(zhí)行更為復(fù)雜和精確的二維或三維數(shù)值模擬，得到更精確的研究成果。

在地質(zhì)構(gòu)造方面，臧紹先等[26]通過有限元法分析俯沖板塊不同參數(shù)變化對俯沖帶負(fù)浮力及其變化過程的影響，揭示了俯沖過程中負(fù)浮力的動態(tài)變化特性。石耀霖等[27]針對活動海嶺俯沖進(jìn)行數(shù)值模擬研究，利用有限元法對該地區(qū)進(jìn)行了熱模擬，結(jié)果表明島弧一定深度下摩擦剪切生熱會形成地溫反轉(zhuǎn)，該過程對島弧火山活動具有重要影響。劉亞靜等[28]采用二維黏彈性數(shù)值模擬技術(shù)，深入研究了俯沖帶深部的應(yīng)力場特性，試驗結(jié)果表明采用適當(dāng)?shù)母_板塊黏度參數(shù)時，橄欖石和尖晶石相變界面下方低黏度區(qū)兩側(cè)會出現(xiàn)應(yīng)力集中區(qū)，相變過渡區(qū)等效黏度對應(yīng)力場的影響不大。唐湘蓉等[29]采用連續(xù)介質(zhì)有限元數(shù)值模擬法，運(yùn)用彈性理論和巖石力學(xué)性質(zhì)，研究裂縫發(fā)育帶的宏觀平面分布，為裂縫預(yù)測提供了新的研究思路。沈海超等[30]利用有限元分析方法討論了斷層如何影響地應(yīng)力場，分析過程中需要考慮多個因素，例如模型需要設(shè)置合適的斷層力學(xué)特性參數(shù)（彈性模量、泊松比等）、斷層介質(zhì)的力學(xué)屬性參數(shù)、斷層的幾何形狀設(shè)置以及邊界條件設(shè)置；研究結(jié)果表明，這些參數(shù)均在不同程度上影響斷層附近的地應(yīng)力場分布，同時導(dǎo)致地應(yīng)力大小方向變化。白玉柱等[31]模擬逆沖斷層在受到垂直于斷層走向的應(yīng)力作用下發(fā)生的逆沖構(gòu)造運(yùn)動，通過使用摩擦接觸單元，并在模型中將斷層面設(shè)計成曲面的幾何形態(tài)，更精確地模擬斷層錯動導(dǎo)致的應(yīng)力分量場改變以及地表在垂直走向和垂直地表上的變化。Davies等[32]提出了一種自適應(yīng)有限元法，改進(jìn)地球動力學(xué)中以對流為主的俯沖帶問題的模擬質(zhì)量。這種自適應(yīng)有限元法在模擬過程中可以根據(jù)誤差自動調(diào)整網(wǎng)格（圖1），提高研究領(lǐng)域流動特征的分辨率。這種自適應(yīng)網(wǎng)格變化由自動非結(jié)構(gòu)化的網(wǎng)格生成器、有限元流求解器以及基于插值的局部誤差估計器實現(xiàn)。Davies等在實驗中驗證了該方法可以有效提高洋脊和俯沖帶模擬的精確度（圖2）。戴黎明等[33]通過有限元法討論亞洲大陸構(gòu)造變形的主要影響因素，根據(jù)當(dāng)今亞洲大陸的構(gòu)造現(xiàn)狀搭建二維有限元模型，模擬主要活動塊體當(dāng)前的構(gòu)造應(yīng)力分布；將模擬結(jié)果與GPS數(shù)據(jù)、震源機(jī)制解和地質(zhì)調(diào)查數(shù)據(jù)進(jìn)行比較，得到主要活動地塊應(yīng)力的分布特征，并對影響構(gòu)造變形的主要因素進(jìn)行了討論分析。袁杰等[34]模擬了斷層自發(fā)破裂的動力學(xué)過程，對滑移弱化摩擦關(guān)系進(jìn)行了改進(jìn)，通過動態(tài)數(shù)值模擬研究了斷層的自發(fā)破裂行為；模擬結(jié)果表明，相較于經(jīng)典的滑移弱化關(guān)系，經(jīng)過優(yōu)化的摩擦關(guān)系可以形成脈沖型（pulse-like）的破裂模式。該實驗還分析了地震破裂過程中的應(yīng)力和應(yīng)變分布，為揭示大地震的破裂機(jī)制提供了新的視角。

在地震研究方面，有限元法主要用于地震期間應(yīng)力場的模擬研究，通過分析地震的相關(guān)參數(shù)和結(jié)果探究地震的成因以及未來預(yù)測。王仁等[35-36]通過數(shù)值模擬技術(shù)對華北地區(qū)地震遷移規(guī)律進(jìn)行了研究，將華北地區(qū)看作多條由主要斷裂帶組成的地質(zhì)構(gòu)造骨架，并在特定深處將模型設(shè)置成理想塑性體，通過分析模擬結(jié)果對華北地區(qū)未來出現(xiàn)強(qiáng)震的可能性作了初步分析；通過模擬實驗，計算出在均勻邊界外力作用影響下華北地區(qū)發(fā)生地震危險的地帶，以及近十二年來歷次大地震后應(yīng)力場的變化，并以此總結(jié)出了一套反演方法。張郢珍等[37]通過有限元法模擬魯南地區(qū)大震的震源應(yīng)力場，研究地震對區(qū)域應(yīng)力場的影響，結(jié)合數(shù)值模擬結(jié)果與實際情況對比，分析震后應(yīng)力應(yīng)變的積累情況，進(jìn)而對魯南地區(qū)未來地震活動的危險性和可能性進(jìn)行了評估和預(yù)測。梅世蓉等[38]利用有限元法對唐山地震的孕震過程進(jìn)行研究，采用二維Maxwell黏彈性模型分析唐山地區(qū)地震的前兆異常機(jī)制和力學(xué)特征，結(jié)果表明唐山地區(qū)可能由于定常邊界力作用，再加上地殼介質(zhì)的非均勻性，導(dǎo)致彈性應(yīng)變得以在局部地區(qū)長期積累形成震源體。王繼存等[39]模擬了川滇菱形塊的構(gòu)造應(yīng)力場，根據(jù)實驗結(jié)果分析討論了該區(qū)域的構(gòu)造應(yīng)力狀況與強(qiáng)震分布之間的聯(lián)系，模擬結(jié)果與實測地震數(shù)據(jù)基本一致，揭示了該地區(qū)構(gòu)造應(yīng)力分布狀況，為強(qiáng)震發(fā)生的可能性分布提供了參考。梁海華等[40]根據(jù)地震震源機(jī)制解和殼幔變形理論建立有限元模型，分析了中國境內(nèi)活動斷裂帶上大震復(fù)發(fā)周期的差異性，結(jié)合古地震資料進(jìn)行區(qū)域板塊邊界力反演，以及構(gòu)造應(yīng)力產(chǎn)生的彈性應(yīng)變能分布特征分析；通過分析應(yīng)變能與板塊邊界距離的衰減關(guān)系，計算得出達(dá)到相同應(yīng)變能所需的累積時間，與已知的斷層帶大震復(fù)發(fā)周期進(jìn)行比較后，得出應(yīng)力衰減的圖像可以有效地解釋不同地區(qū)大震復(fù)發(fā)周期的顯著差異。Parsons[41]利用三維有限元模擬技術(shù)研究了1906年舊金山大地震后圣安地列斯斷層系統(tǒng)的應(yīng)力恢復(fù)問題。相較于前人，Parsons的模型僅關(guān)注斷層的構(gòu)造應(yīng)力，顯示了更為均勻的地殼應(yīng)力分布。研究結(jié)果為理解大地震后應(yīng)力陰影對后續(xù)地震活動的影響提供了新的模擬技術(shù)支持，同時提供更為準(zhǔn)確的數(shù)據(jù)參考，對舊金山灣區(qū)的未來地震的危險性和可能性進(jìn)行了評估。朱守彪等[42]采用了有限單元法模擬蘇門答臘俯沖帶大地震的發(fā)生過程，通過建立兩種描述了影響速率中黏著和滑移摩擦接觸運(yùn)動狀態(tài)的數(shù)學(xué)模型，采用時間積分方法對模型求解來保證力學(xué)變化的穩(wěn)定性，模擬結(jié)果表明俯沖帶上大地震的發(fā)生必須滿足大規(guī)模、均勻介質(zhì)和相同摩擦系數(shù)的區(qū)域這3個條件。李平恩等[43]利用三維黏彈性有限元模型技術(shù)，依據(jù)青藏高原巴顏喀拉塊體及其周緣地區(qū)的地質(zhì)構(gòu)造特征搭建深部速度模型，模型設(shè)置考慮了巴顏喀拉塊體區(qū)域地質(zhì)構(gòu)造、活動塊體、活動斷裂帶以及邊界斷裂帶的分層。Qu等[44]基于庫侖破裂應(yīng)力，利用數(shù)值模擬建立模型，討論了應(yīng)力變化與強(qiáng)震之間的關(guān)系、強(qiáng)震之間的互相作用，以及長時間的構(gòu)造加載如何影響強(qiáng)震的發(fā)生；模擬結(jié)果表明2011年日本Mw 9.0東北地震和2008年汶川Mw 7.9地震對華北地區(qū)的地殼構(gòu)造活動產(chǎn)生了顯著影響。Bahuguna等[45]研究了印度板塊的板內(nèi)應(yīng)力分布，實驗采用三維力學(xué)模型，設(shè)置模型過程中充分考慮了印度次大陸和洋區(qū)19個地質(zhì)區(qū)域的彈性特性，建模采用ABAQUS軟件，最終的模擬結(jié)果與板塊速度一致。Bahuguna等的研究能夠提升地震和地質(zhì)研究的解釋性，建模方法可用于識別地震活躍區(qū)域并進(jìn)一步對地震危險性進(jìn)行評估，通過解釋應(yīng)變率、變形率和應(yīng)力分布數(shù)據(jù)，該研究為印度板塊地震活躍區(qū)域的活動板內(nèi)變形提供了新的解釋。Karabulut等[46]通過有限元法對Main Marmara Fault（MMF）進(jìn)行了三維構(gòu)造建模，并評估不同網(wǎng)格類型在揭示地震潛力方面的準(zhǔn)確性。三維模型數(shù)據(jù)來源于位于MMF周圍GNSS站點的速度數(shù)據(jù)，并選擇ANSYS有限元軟件進(jìn)行模擬，檢驗使用不同網(wǎng)格類型構(gòu)建模型的準(zhǔn)確性。通過比較模型中斷層滑移虧缺值的最低和最高均方根誤差，得出了可能的地震矩震級差異在0.010～0.014之間。該研究側(cè)重點在于進(jìn)行有限元分析時選擇合適網(wǎng)格類型的重要性，同時實驗結(jié)果也為地震潛力評估和風(fēng)險評估提供了新方法。馬峰等[47]利用數(shù)值模擬方法對雄安新區(qū)容城地?zé)崽锼E縣系碳酸鹽巖的熱儲進(jìn)行模擬，通過實驗驗證在不同的采灌條件下熱儲的溫度和壓力反應(yīng)。Xu等[48]提出了一種多尺度擴(kuò)展有限元法（MS-XFEM），適用于模擬地質(zhì)構(gòu)造中多重斷裂的傳播。該方法通過局部更新基函數(shù)來追蹤斷裂尖端的動態(tài)變化，無須額外的自由度，因此有效降低了計算成本。MS-XFEM通過結(jié)合預(yù)處理的GMRES（generalized minimum residual）迭代策略控制最終解的精度，確保斷裂路徑的正確預(yù)測。通過測試結(jié)果來看，MS-XFEM方法能夠有效地模擬復(fù)雜的斷裂模式，即使在較大的誤差容限下也能保持較高的預(yù)測準(zhǔn)確性，為斷裂模擬研究提供了一種新的計算方法。Chen等[49]提出了一種在交錯斷層帶中生成層狀地質(zhì)體三維有限元網(wǎng)格的方法。在三維地質(zhì)建模過程中，斷層帶空間關(guān)系復(fù)雜，需要模型具備一定的精確性。該研究通過分析交錯斷層帶中不同網(wǎng)格類型之間的空間關(guān)系，提出了一種封閉流形處理方法，并通過該方法建立交錯斷層帶的封閉流形空間網(wǎng)格模型來劃分地質(zhì)體3D網(wǎng)格。

常見有限單元法對比見表1。

3.2 有限差分法

有限差分法是最早的數(shù)值模擬算法，具有簡單直觀的優(yōu)勢，其基本思想就是用泰勒級數(shù)展開式在覆蓋計算域的網(wǎng)格或節(jié)點上，有限差分替換無限微分，利用離散數(shù)值集合來近似原有的連續(xù)函數(shù)場。一般情況下需要對邊界條件作特殊處理。該方法在處理規(guī)則幾何體時優(yōu)勢明顯。

有限差分法最早由歐拉（L. Euler）[50]在1768年提出，用于求解一維問題的差分格式。1908年，龍格（C. Runge）[51]將差分法延展到了二維的問題上。Alterman等[52]于1968年提出使用差分算子代替波動方程中的微分項，將波動方程離散化，并成功實現(xiàn)了2D層狀模型中的彈性波模擬。這一開創(chuàng)性的工作采用顯式有限差分法得到了二階彈性波方程對應(yīng)的離散數(shù)值解。隨后學(xué)者們對該方法的原理和精度進(jìn)行了深入分析和研究，并進(jìn)一步將其推廣應(yīng)用到聲波介質(zhì)和非均勻介質(zhì)中。1986年，Virieux[53]基于交錯網(wǎng)格進(jìn)行了非均勻介質(zhì)中P-SV（pressure-shear）波有限差分模擬，指出基于完全交錯網(wǎng)格可將彈性波模擬的方法推廣到S（shear）波速度為零的液體介質(zhì)，且不需要對液-固界面做特殊處理，并通過解析解對比以及相應(yīng)的數(shù)值模型驗證了此結(jié)論；Dablain[54]指出可以使用二階中心差分格式的高階來近似提高模擬精度和效率；Bayliss等[55]推導(dǎo)并實現(xiàn)了空間四階、時間二階差分格式下的彈性波正演模擬，同時將四階格式與二階格式的精度進(jìn)行了對比分析；Levander[56]基于Virieux給出的交錯網(wǎng)格推導(dǎo)了精確空間四階、時間二階差分算子以及二維P-SV波有限差分格式，并指出每個波長內(nèi)至少需要5個網(wǎng)格才能滿足精度需求；Jastram等[57]提出了一種垂向變間距的網(wǎng)格格式，增加了靈活性并提高了計算效率；在此基礎(chǔ)上，F(xiàn)alk等[58]又提出一種變時間步長的有限差分格式，大大降低了計算成本；Graves[59]提出了一種交錯網(wǎng)格有限差分技術(shù)，并應(yīng)用于三維彈性介質(zhì)模型中的地震波傳播，提出了具有有效材料參數(shù)的三維四階速度-應(yīng)力格式，同時介紹了一種內(nèi)存優(yōu)化算法，使得在標(biāo)準(zhǔn)工作站上進(jìn)行大規(guī)模三維模擬成為可能。

20世紀(jì)60年代后期，有限差分法才被應(yīng)用到地震波傳播的模擬當(dāng)中，這也是有限差分法最先在地學(xué)中應(yīng)用的開端。1972年，Boore[60]提出了一種有限差分法，并用于模擬在非均質(zhì)介質(zhì)中的P-SV波傳播。該方法在原有SH（shear horizontal）波傳播模型基礎(chǔ)上，通過在離散網(wǎng)格中使用速度和應(yīng)力進(jìn)行擴(kuò)展，驗證了穩(wěn)定性條件與泊松比無關(guān)，同時也表明了相同的代碼可以用于彈性介質(zhì)和液體介質(zhì)，其中S波速度為零，且液體-固體界面無須做特殊處理。Kelly等[61]研究如何利用有限差分法生成合成地震圖，利用該方法在復(fù)雜的地下結(jié)構(gòu)幾何和任意源-接收器距離條件下產(chǎn)生合成地震圖。該研究不僅提供了對地震圖解釋的見解，同時也對利用有限差分法模擬過程中各種現(xiàn)象做出了論述，如網(wǎng)格色散、模型邊緣的人工反射以及空間和時間采樣間隔的選擇。Madariaga[62]于1976年研究了通過數(shù)值模擬方法模擬平面圓形摩擦斷層模型的動力學(xué)，詳細(xì)討論了遠(yuǎn)場輻射，并區(qū)分了三個頻譜區(qū)域：通常的恒定低頻水平、受斷層尺寸控制的中間區(qū)域，以及高頻漸近線。由此提出了拐角頻率與遠(yuǎn)場位移脈沖寬度成反比的觀點，這為斷層尺寸的測量提供了物理依據(jù)，斷層體大小與前者存在定量關(guān)系。高孟潭等[63]利用三維有限差分法對北京地區(qū)地震動進(jìn)行了模擬研究，并通過模擬實驗討論了盆地效應(yīng)對地震動的影響。Virieux[53]提出了一種用于模擬非均質(zhì)介質(zhì)中P-SV波傳播的有限差分法，該方法解決了P-SV波傳播的源建模和表面波（Rayleigh波）建模兩個問題。并將該方法的模擬結(jié)果與簡單解析解進(jìn)行了比較，驗證了方法的有效性，最后展示了在不同泊松比下液體-固體界面的復(fù)雜幾何形狀的模擬結(jié)果。Hernlund等[64]基于三維地幔對流模型提出了求解poloidal速度勢的有限差分法。Igel等[65]于1992年首次推導(dǎo)出各向異性介質(zhì)中交錯網(wǎng)格有限差分形式，并在文章中詳細(xì)討論了該形式模擬精度及影響因素，這種有限差分法可以實現(xiàn)強(qiáng)各向異性介質(zhì)中的地震波模擬。Pitarka[66]將非均勻間距的網(wǎng)格應(yīng)用到三維彈性介質(zhì)中，大大降低了內(nèi)存消耗，同時提高了三維數(shù)據(jù)的處理效率。Bohlen等[67]也提出了一種旋轉(zhuǎn)交錯網(wǎng)格（RSG）有限差分技術(shù)，并將這種方法成功應(yīng)用到各向異性黏彈性介質(zhì)的地震波模擬中，最終得到了較為準(zhǔn)確的模擬結(jié)果。

進(jìn)入21世紀(jì)，有限差分法在模擬近地表復(fù)雜介質(zhì)方面的應(yīng)用得到了進(jìn)一步的發(fā)展，研究者開始關(guān)注如何在有限差分法框架下處理復(fù)雜地表地形的起伏問題。近年來，有限差分法在地學(xué)領(lǐng)域的應(yīng)用繼續(xù)延展，當(dāng)前已廣泛應(yīng)用于地震波傳播、地下結(jié)構(gòu)分析、資源勘探等方面。李建平等[68]采用球坐標(biāo)系下的交錯網(wǎng)格有限差分法進(jìn)行地磁測深的三維正演模擬。該方法主要基于Maxwell方程的積分形式，數(shù)學(xué)模型離散后的方程組采用PARDISO求解器處理，避免迭代求解中的散度校正問題。該研究中還與已有的有限元法和有限差分法進(jìn)行了比較，該方法在一維層狀模型和雙半球模型中顯示出較高的計算精度，相對誤差優(yōu)于5%。韓麗等[69]利用基于Maxwell全方程的頻域有限差分法對地下多孔介質(zhì)中震電效應(yīng)進(jìn)行數(shù)值模擬，通過精確模擬地震波與電磁波場的相互作用，研究不同介質(zhì)分界面處的電磁波場反射和透射現(xiàn)象，準(zhǔn)確識別異常體位置，有助于更好地理解震電效應(yīng)產(chǎn)生的機(jī)制。張志佳等[70]利用有限差分法和起伏多重變加密網(wǎng)格技術(shù)對復(fù)雜地表波動方程進(jìn)行數(shù)值模擬，誤差控制在10-12內(nèi)，計算效率較高，可以很好地適應(yīng)實際復(fù)雜地質(zhì)條件。Gerya等[71]提出了一種新的交錯網(wǎng)格有限差分法，對網(wǎng)格分辨率轉(zhuǎn)換過程中內(nèi)部的速度節(jié)點施加了應(yīng)力守恒約束，使網(wǎng)格在離散化的過程中得以保持矩陣稀疏性，大大降低了計算成本，提高了求解效率。該方法適用于強(qiáng)烈可變黏度情況的模擬，Gerya等利用該方法進(jìn)行了巖石圈結(jié)構(gòu)模擬（圖3）和行星尺度模擬（圖4），結(jié)果顯示在黏度發(fā)生顯著變化的區(qū)域或復(fù)雜計算域地區(qū)，該方法會自動增加網(wǎng)格密度，保證計算效率同時也保證了準(zhǔn)確率。目前研究者們?nèi)栽诓粩喔倪M(jìn)算法，以提高計算效率和精度，將有限差分法和其他數(shù)值方法結(jié)合，如有限元法和邊界元法等，以此來解決更加復(fù)雜的地學(xué)問題。Chávez-Negrete等[72]使用廣義有限差分法模擬地下水流，提出了一種新的界面公式來應(yīng)對非飽和多孔層狀介質(zhì)中的地下水流現(xiàn)象，確保了在不同介質(zhì)層之間的界面節(jié)點處水的平衡。

常見有限差分法對比見表2。

3.3 邊界元法

邊界元法是一種求解邊界積分方程中邊界值或初值問題的數(shù)值方法框架，主要用于求解邊界條件對解起主導(dǎo)性作用的物理問題。邊界元法將問題域的求解過程轉(zhuǎn)換為邊界上的積分方程，并應(yīng)用數(shù)值離散化手段，將連續(xù)的問題轉(zhuǎn)化為線性代數(shù)方程組的形式。與其他基于問題域的數(shù)值方法相比較，邊界元法相當(dāng)于將問題降維，更容易生成網(wǎng)格。邊界元法擅長求解大規(guī)模邊界問題，例如地殼和地幔建模、計算流體動力學(xué)和熱傳導(dǎo)問題等。

邊界元法的理論根基在于積分方程的數(shù)學(xué)理論，其中基本解（Green函數(shù)）為從邊界到問題域內(nèi)部任意點的影響關(guān)系提供了數(shù)學(xué)描述。邊界元法能夠有效地處理無限或半無限域問題，該問題在傳統(tǒng)的基于域的數(shù)值方法中需要在無限處引入額外的條件，處理通常較為復(fù)雜。此外，邊界元法在處理邊界條件變化或邊界形狀不規(guī)則的問題時具有很好的靈活性和高效性。但同時也需要注意邊界元法的精確性很大程度上依賴于邊界離散化的質(zhì)量，以及數(shù)值積分方案是否嚴(yán)格。

邊界元法的出現(xiàn)相較于有限元法晚了約20年，最早是Jaswon[73]于1963年討論了Fredholm積分方程與Green邊界公式在位勢理論和彈性力學(xué)的應(yīng)用。同時該方法的提出也為困難邊界問題，如非彈性動力學(xué)問題以及其他梁或板組成的結(jié)構(gòu)提供了解決方法。Rizzo[74]提出一種基于矢量邊界公式的方法，無須復(fù)雜映射函數(shù)就能處理多連通域問題。這種方法展示了邊界元法處理彈性靜力學(xué)問題的能力。Sweedlow等[75]利用邊界元法對三維彈性標(biāo)本進(jìn)行了數(shù)值應(yīng)力分析，展示了裂紋尖端附近應(yīng)力分布及其三維特性。同年，Cruse等[76]又進(jìn)一步探討了邊界元法在三維彈塑性流動問題中的應(yīng)用，通過擴(kuò)展Somigliana恒等式，在邊界上評估材料關(guān)系降低維度，這種方法適用于一般性彈塑性流動理論。Mendelson[77]研究了邊界元在彈塑性扭轉(zhuǎn)問題的應(yīng)用，并與有限差分法進(jìn)行了比較，精度與之類似，收斂速度更快。Yin等[78]提出一種新的正則化公式，有效提高了數(shù)值解的準(zhǔn)確性和穩(wěn)定性，其通過直接法和間接法將二維流固耦合問題轉(zhuǎn)化為三組不同的邊界積分方程，并引入邊界算子進(jìn)行求解。其中直接法基于Green公式和基本解來表示經(jīng)典解，間接法基于勢層理論，這兩種方法都可以將問題轉(zhuǎn)化為邊界值問題。數(shù)值模擬結(jié)果表明該方法可以有效解決流固耦合問題。Quevedo等[79]利用快速多極邊界元法模擬了巖石圈-地表相互作用動力演化形成的洋中脊地形演變（圖5）。該方法融合了邊界元和多極方法能夠提高計算應(yīng)力和應(yīng)變的效率，減少了所需的網(wǎng)格元素數(shù)量，使求解問題所需的內(nèi)存和時間需求大大降低。該方法適用于處理地殼斷層和地殼內(nèi)物質(zhì)轉(zhuǎn)換的復(fù)雜性問題，能夠處理大規(guī)模的地球動力學(xué)模型（如全球地殼動力學(xué)模擬）。

我國關(guān)于邊界元的研究始于20世紀(jì)70年代后期，馮康教授在1978年根據(jù)橢圓邊值的特性，首次提出邊界歸化思想；同一時期清華大學(xué)杜慶華教授、姚振漢教授，以及嵇醒、申光憲等學(xué)者也開始從事邊界元的研究，并在彈性分析等領(lǐng)域作出重要貢獻(xiàn)[80]。在20世紀(jì)的80年代，張永元、祝家麟、黃玉盈、胡海昌、王有成等為代表的科學(xué)家對邊界元法在流固相互作用、彈塑性接觸以及板殼結(jié)構(gòu)分析等領(lǐng)域進(jìn)行了廣泛應(yīng)用和研究[81]。

在20世紀(jì)90年代，邊界元法的研究繼續(xù)深化，被應(yīng)用于更復(fù)雜的領(lǐng)域，如復(fù)雜彈性力學(xué)、多相流、復(fù)合材料結(jié)構(gòu)分析、非線性彈塑性和流固耦合等，算法精度也得到進(jìn)一步提升[82]。

進(jìn)入21世紀(jì)后，邊界元法在地球動力學(xué)的非線性問題求解中得到了廣泛的采納和應(yīng)用。金朝嵩[83]通過Sobole空間框架對一般算子方程的Galerkin逼近進(jìn)行了深入分析，并給出了后驗誤差估計的結(jié)果。針對以有限平面為屏蔽物的聲散射問題，該研究為自適應(yīng)邊界元法的解法提供了具體的后驗誤差估計表達(dá)式。該研究能夠有效地處理邊界積分方程中的奇點問題，提高數(shù)值解的精度和效率。2003年，吳云等[84]提出了一種基于邊界元的非連續(xù)（塊體系統(tǒng)）形變反分析法，采用邊界元法和最小二乘原理，以邊界單元的位移和表面力作為未知參數(shù)，以觀測到的位移作為觀測值，來確定邊界上的相對運(yùn)動、位移和應(yīng)力場。該研究討論了如何確定實驗建模過程中為保證系統(tǒng)整體的正確性和合理性所需的應(yīng)力參數(shù)和其他約束條件，以及如何判斷塊體間接觸狀態(tài)。2010年，焦龍等[85]實現(xiàn)了對川西地區(qū)形變場與應(yīng)變場的研究，分析了巴顏喀拉地塊和四川地塊在汶川地震前的地殼形變特征。本次研究將GPS位移數(shù)據(jù)作為輸入，計算塊狀邊界應(yīng)力數(shù)值，以及內(nèi)部應(yīng)力應(yīng)變和位移。研究結(jié)果表明研究地塊南北向的張應(yīng)變作用以及東西向應(yīng)變能的積累可能是觸發(fā)汶川地震的關(guān)鍵因素。

常見邊界元法對比見表3。

3.4 有限體積法

有限體積法最早可以追溯到20世紀(jì)50年代，該方法的主要應(yīng)用領(lǐng)域為航空航天領(lǐng)域的計算流體力學(xué)。

有限體積法的基本思想與有限差分類似，均以問題的積分公式入手，將計算域轉(zhuǎn)化為一系列控制體積，這些控制體積通常與網(wǎng)格單元相關(guān)聯(lián)；之后對每個控制體積應(yīng)用守恒定律進(jìn)行離散化得到一組離散方程。這種方法使得離散方程不僅具有明確的物理意義，同時也確保了物理量的守恒特性。

有限體積法在地學(xué)領(lǐng)域的應(yīng)用雖然起步較晚，但發(fā)展迅速。首次被引入地學(xué)領(lǐng)域是為了解決地下水流動和污染物運(yùn)移的問題。隨著時間的推移，有限體積法在地學(xué)中的應(yīng)用被廣泛擴(kuò)展到多個領(lǐng)域，包括地質(zhì)建模、地震模擬、滑坡分析、海洋流體動力學(xué)和氣候變化研究等。

最早使用有限體積法離散化的計算代碼為Ogawa 等[86]在1991年開發(fā)的三維笛卡爾模型，該模型采用了迭代求解器。之后，Tackley[87]利用三維笛卡爾代碼Stag3D實現(xiàn)了該網(wǎng)格上的速度-壓力多網(wǎng)格求解器，并將該方法應(yīng)用于模擬可壓縮對流，目前該代碼已在眾多科學(xué)研究中得到推廣應(yīng)用。Stag3D軟件采用了MPDATA平流方法（一種有限體積法）。

隨后開發(fā)的若干計算軟件或者代碼雖然同樣采用了基礎(chǔ)的有限體積/有限差分多網(wǎng)格方法，但都根據(jù)具體研究問題進(jìn)行了相應(yīng)的改進(jìn)。Gerya等[88]開發(fā)的二維代碼采用了相同的錯開網(wǎng)格，但包含了許多額外的復(fù)雜物理特性，如彈性、自由表面和多種成分物質(zhì)的追蹤等。Choblet等[89]采用非多網(wǎng)格迭代求解器實現(xiàn)了有限體積法在經(jīng)度-緯度-半徑網(wǎng)格上三維球面幾何的應(yīng)用。使用立方球體網(wǎng)格、陰陽網(wǎng)格以及螺旋網(wǎng)格實現(xiàn)了多網(wǎng)格版本。

Aziz等[90]在石油儲層模擬中應(yīng)用有限體積法，介紹了黑油模型和組分模型的數(shù)值求解，為油氣資源的高效開發(fā)提供了重要工具。Jung等[91]采用有限元和有限體積混合方法，使用隱式壓力-顯式飽和度數(shù)值方案處理了流體相之間的毛細(xì)力和其他非線性效應(yīng)，模擬了石油在盆地中的長距離遷移和積累過程（圖6）。該方法能有效地處理地質(zhì)結(jié)果中非均質(zhì)性和流體動力學(xué)的復(fù)雜性，為油氣勘探提供了有效模擬方法。馬瑞杰等[92]提出了一種新的有限體積法數(shù)學(xué)模型。該模型根據(jù)水流連續(xù)性原理、質(zhì)量守恒定律和達(dá)西定律，其中巖溶裂隙介質(zhì)被定義為雙重孔隙度介質(zhì)，以此建立非穩(wěn)定流雙重介質(zhì)數(shù)學(xué)模型對裂隙地下水流動問題進(jìn)行求解。Jenny等[93]開發(fā)一種高效的有限體積法用于多相流和傳輸?shù)亩喑叨饶M，處理非線性流動和傳輸問題。

利用有限體積法對復(fù)雜地質(zhì)構(gòu)造中的地下水流動和污染物運(yùn)移進(jìn)行了數(shù)值模擬實驗，為環(huán)境保護(hù)和管理提供了新的技術(shù)手段和研究視角。Cui等[95]采用有限體積法對地?zé)醿又械乃疅狁詈线^程進(jìn)行了模擬，模擬結(jié)果展示了該方法在地?zé)豳Y源評估和開發(fā)中的應(yīng)用。Huang等[96]對裂隙巖中的地下水流動和溶質(zhì)運(yùn)移進(jìn)行了數(shù)值模擬，結(jié)果表明有限體積法能高效準(zhǔn)確地對裂隙介質(zhì)中的水流和溶質(zhì)遷移進(jìn)行研究。

3.5 無網(wǎng)格法

無網(wǎng)格法克服了傳統(tǒng)方法對網(wǎng)格的依賴性，通過一系列散布在求解域內(nèi)的節(jié)點來近似表示未知函數(shù)。該方法適用于涉及網(wǎng)格不規(guī)則、邊界移動、形變較大及高維復(fù)雜幾何域的問題。這種特性使得無網(wǎng)格法在地學(xué)領(lǐng)域的復(fù)雜問題建模中具有較高的應(yīng)用價值。

無網(wǎng)格數(shù)值技術(shù)的研究最早出現(xiàn)于20世紀(jì)70年代，當(dāng)時主要圍繞任意網(wǎng)格有限差分法進(jìn)行研究。1977年Lucy[97]提出了一種創(chuàng)新的粒子模擬方法——光滑粒子流體動力學(xué)法，該方法被應(yīng)用于天體物理學(xué)中的復(fù)雜現(xiàn)象，如行星碰撞和星系演化。隨著研究的深入，光滑粒子流體動力學(xué)法的應(yīng)用逐漸擴(kuò)展到流體動力學(xué)等更廣泛的領(lǐng)域。盡管光滑粒子流體動力學(xué)法最初在處理邊界條件和非均勻粒子分布方面存在一定的局限性，但后續(xù)的研究工作對其不斷進(jìn)行了有效的改進(jìn)和優(yōu)化。

在無網(wǎng)格法的發(fā)展過程中，Nayroles等[98]于1992年提出了離散單元法，該方法通過加權(quán)最小二乘多項式擬合實現(xiàn)了局部近似解的構(gòu)建。Belytschko 等[99]在1994年提出了無單元Galerkin法，該方法采用拉格朗日乘子技術(shù)確保滿足強(qiáng)邊界條件，顯著提升了數(shù)值解的收斂速度。無單元Galerkin法在邊界條件處理方面表現(xiàn)優(yōu)越的同時，其較高的計算成本和對背景網(wǎng)格的依賴性成為這種方法被廣泛應(yīng)用的制約因素。

Liu等[100]在1995年提出了再生核粒子法，該方法通過修正函數(shù)保證了邊界條件的滿足，又進(jìn)一步提升了數(shù)值解的精度。與此同時提出有限點法和徑向基點插值法等技術(shù)，它們由于無須背景網(wǎng)格的特性在流體力學(xué)和跨學(xué)科數(shù)值模擬中得到了廣泛應(yīng)用。

徑向基函數(shù)法的研究最早來源于Hardy[101]在1971年提出的無網(wǎng)格法，通過稀疏數(shù)據(jù)描述復(fù)雜的地形表面。這種方法的優(yōu)勢在于構(gòu)建形函數(shù)無須考慮復(fù)雜情況，只需要考慮離散點間的距離，有助于高維問題的降維處理。隨后，引入正定緊支徑向函數(shù)和緊支徑向基函數(shù)，有效解決了無網(wǎng)格法中系數(shù)矩陣的稀疏性問題，也為后續(xù)大型問題的求解奠定了基礎(chǔ)。局部微分求積的徑向基函數(shù)法和結(jié)合高斯函數(shù)的方法能夠進(jìn)一步增強(qiáng)數(shù)值解的穩(wěn)定性和精度。

無網(wǎng)格數(shù)值技術(shù)已經(jīng)形成了一個包含光滑粒子流體動力學(xué)法、離散元法、無單元Galerkin法、再生核粒子法、有限點法、徑向基點插值法等多種方法的豐富體系，這些技術(shù)在各自的應(yīng)用領(lǐng)域內(nèi)展現(xiàn)出獨特的優(yōu)勢和潛力。盡管無網(wǎng)格技術(shù)在理論完備性、計算效率和實用性方面仍需進(jìn)一步的研究和改進(jìn)，但在計算力學(xué)領(lǐng)域的研究熱潮和應(yīng)用前景不容忽視。學(xué)術(shù)界對此領(lǐng)域進(jìn)行了深入研究，如張雄等[102]的著作和Nguyen等[103]的綜述文章，為無網(wǎng)格技術(shù)的理解和應(yīng)用提供了寶貴的知識和資源。

焦玉玲等[104]利用無網(wǎng)格法實現(xiàn)了地下水水位的預(yù)測，分析了鞍山市首山區(qū)水文地質(zhì)條件對水位變化的影響。該研究采用移動最小二乘理論的無網(wǎng)格法，通過構(gòu)建場函數(shù)、基函數(shù)和權(quán)函數(shù)的選取，形函數(shù)及其導(dǎo)數(shù)的計算，建立了雙層滲流二維平面系統(tǒng)的數(shù)學(xué)模型。模擬預(yù)測得到的結(jié)果與實際水位變化規(guī)律基本一致，驗證了無網(wǎng)格法在水量模擬計算中的有效性。王月英[105]采用無單元Galerkin法進(jìn)行了地震波正演模擬，并對其影響因素進(jìn)行了深入分析。該方法的優(yōu)勢在于通過局部支撐域上的權(quán)函數(shù)實現(xiàn)最佳逼近，突破了傳統(tǒng)網(wǎng)格剖分的局限。但需要注意的是基函數(shù)的選擇對模擬精度和計算負(fù)荷有顯著影響，基函數(shù)階數(shù)越高，精度越高，導(dǎo)致運(yùn)算量越大；同時節(jié)點的分布和編號方式會影響數(shù)據(jù)量和運(yùn)算效率。該研究展示了無單元Galerkin法在地震波模擬中應(yīng)用的潛力。李坤等[106]采用無網(wǎng)格法求解二階時域波動方程，依據(jù)徑向基函數(shù)逼近空間導(dǎo)數(shù)以及Crank-Nicolson時間離散化技術(shù)進(jìn)行數(shù)值模擬。該研究采用的方法可提供與有限元相似的精度，同時與有限元法相比采用了更少的數(shù)據(jù)點。該研究證明了無網(wǎng)格法在處理復(fù)雜幾何問題時的優(yōu)勢，相較于其他數(shù)值計算方法簡化了計算流程，避免了網(wǎng)格劃分和重構(gòu)的步驟。由于該研究方法采用的徑向基函數(shù)與問題維度無關(guān)，所以可擴(kuò)展到更高維度問題的求解。Wang等[107]提出了一種基于徑向基函數(shù)的無網(wǎng)格法用于求解逆波傳播問題。該方法結(jié)合了配點法和徑向基函數(shù)，通過Tikhonov正則化技術(shù)和L-曲線準(zhǔn)則處理噪聲測量數(shù)據(jù)，提高了解決問題的穩(wěn)定性。研究對比分析了傳統(tǒng)有限元法和流行的弱形式無網(wǎng)格法，相較而言該方法具有更少的色散誤差，并通過穩(wěn)定性分析預(yù)測了顯式時間積分的關(guān)鍵時間步長。馮德山等[108]基于滑動最小二乘法擬合場函數(shù)采用無單元Galerkin法對探地雷達(dá)進(jìn)行正演模擬。無單元Galerkin法避免了復(fù)雜的網(wǎng)格剖分，簡化計算流程，保持了高計算精度和解的高次連續(xù)性。該研究對無單元Galerkin法處理強(qiáng)加邊界條件的策略進(jìn)行了創(chuàng)新，采用罰因子法和透射邊界條件提高模擬的準(zhǔn)確性。通過與基于線性插值的有限元法正演剖面圖的對比，證明無單元Galerkin法在相同節(jié)點數(shù)條件下具有更高的精度，能更好地解釋探地雷達(dá)數(shù)據(jù)。無網(wǎng)格法也被用于模擬地震波與建筑物、橋梁等結(jié)構(gòu)的相互作用。例如，嵇艷鞠等[109]將無網(wǎng)格法應(yīng)用于起伏地形各向異性2D大地電磁響應(yīng)模擬，分析了起伏地形對電磁響應(yīng)的影響，為電磁法勘探中的地形校正和數(shù)據(jù)解釋提供了新的研究方法。李俊杰等[110]利用無網(wǎng)格點插值法對大地電磁場響應(yīng)進(jìn)行數(shù)值模擬，研究分析了無網(wǎng)格點插值法在處理復(fù)雜地質(zhì)模型時的優(yōu)勢，特別是在與有限元法和無單元Galerkin法的比較中展示了無網(wǎng)格點插值法在計算精度和效率方面的顯著提升。魏德敏等[111]將無單元Galerkin法應(yīng)用于土質(zhì)邊坡的彈塑性穩(wěn)定分析，通過與有限元法計算結(jié)果的對比分析，發(fā)現(xiàn)使用無單元Galerkin法計算彈塑性問題收斂速度快、穩(wěn)定性高。該研究為彈塑性問題的求解提供了一種新的數(shù)值計算方法。McDougall等[112]利用無網(wǎng)格法提出了一種可以用來動態(tài)分析快速流動滑坡、泥石流和雪崩的數(shù)值模型，該模型基于平滑粒子流體動力學(xué)（SPH）模擬復(fù)雜三維地形上的流動。該方法適用于非靜水壓力和各向異性內(nèi)部應(yīng)力分布，使得快速流動現(xiàn)象的動態(tài)特性能夠在模擬中被展示出來，同時能夠消除長距離位移網(wǎng)格扭曲問題。他們在實驗室水槽實驗中分析了真實巖石雪崩案例（圖7），從而證明了該模型模擬實際滑坡、巖屑流和雪崩問題的能力。Lin等[113]提出了一種改進(jìn)的無網(wǎng)格數(shù)值流形方法。該方法采用移動最小二乘插值，解決了傳統(tǒng)方法在模擬復(fù)雜邊界滲透問題時出現(xiàn)的低精度和網(wǎng)格依賴性問題。Cao等[114]提出了一種無網(wǎng)格數(shù)值流形方法的改進(jìn)版，用于邊坡穩(wěn)定性分析。該研究通過結(jié)合強(qiáng)度折減法和Mohr-Coulomb（M-C）準(zhǔn)則的多屈服面控制方程，解決了有限元法在計算具有軟弱夾層的邊坡時對網(wǎng)格數(shù)量和非線性方程求解算法的依賴性和敏感性問題。研究中引入了主應(yīng)力空間中的子空間跟蹤方法來處理M-C多屈服面的角點問題，并結(jié)合k（k為聚類中心數(shù)）均值聚類算法和小波變換，開發(fā)了一種智能且高效的臨界滑移面自動提取方法。

常見無網(wǎng)格法對比見表4。

3.6 數(shù)值計算方法對比

在地球動力學(xué)數(shù)值模擬研究中，不同的數(shù)值計算方法各有優(yōu)缺點和適用領(lǐng)域。選擇合適的數(shù)值計算方法對于準(zhǔn)確模擬地球內(nèi)部過程至關(guān)重要。例如有限元法處理復(fù)雜幾何形狀和材料屬性具有優(yōu)勢，有限差分法在規(guī)則網(wǎng)格和簡單邊界問題時展現(xiàn)出較高的計算效率，有限體積法更擅長處理流體動力學(xué)，邊界元法能高效處理邊界問題，無網(wǎng)格法則能夠靈活處理大變形和復(fù)雜運(yùn)動。表5對幾種主要的數(shù)值計算方法進(jìn)行了對比分析，以幫助研究人員更好地理解這些方法的特點和適用范圍，從而在實際研究中根據(jù)實際情況做出最優(yōu)選擇。

3.7 數(shù)值模擬軟件、代碼

基于上述各種數(shù)值計算方法，研究者們在地球動力學(xué)數(shù)值模擬領(lǐng)域內(nèi)開發(fā)了眾多軟件和代碼，每個工具都具有獨特的應(yīng)用領(lǐng)域和特點，以滿足不同的研究需求。這些軟件和代碼在地球動力學(xué)研究中扮演著至關(guān)重要的角色，它們使研究人員能夠?qū)⒏嗟木性诳茖W(xué)問題而非算法上，進(jìn)而幫助研究人員模擬和理解復(fù)雜的地球演化過程，為推動地球科學(xué)發(fā)展做出了卓越的貢獻(xiàn)。表6展示常見數(shù)值模擬軟件、代碼的特點和應(yīng)用領(lǐng)域。

4 總結(jié)及展望

本文綜述了數(shù)值模擬在地球動力學(xué)研究中的應(yīng)用現(xiàn)狀與進(jìn)展，分析了各種數(shù)值方法在模擬地球系統(tǒng)動力學(xué)過程中的優(yōu)勢與局限。通過對比有限元法、有限差分法、有限體積法、邊界元法和無網(wǎng)格法等主要數(shù)值計算方法，我們可以發(fā)現(xiàn)每一種方法都有其特定的適用場景和要求。

有限元法因其高度的靈活性和適應(yīng)性，在處理復(fù)雜地質(zhì)結(jié)構(gòu)和邊界條件方面表現(xiàn)出色，尤其適用于大尺度和長時間尺度的模擬；有限差分法和有限體積法在流體流動和熱傳導(dǎo)問題上具有高效率，適合于規(guī)則幾何和均質(zhì)介質(zhì)的模擬；邊界元法在處理無限或半無限域問題時具有明顯優(yōu)勢，能有效減少計算量；無網(wǎng)格法則在處理大變形和破壞問題時顯示出其強(qiáng)大的能力，尤其是在不需要傳統(tǒng)網(wǎng)格的復(fù)雜幾何形狀模擬中。

未來的發(fā)展趨勢將聚焦于以下幾個方面：第一，多物理場耦合模擬技術(shù)的進(jìn)步將推動地球系統(tǒng)模擬的深入發(fā)展，這將有助于更全面地理解地球系統(tǒng)的復(fù)雜性，揭示其內(nèi)在機(jī)理和演變規(guī)律。通過耦合多種物理場，模擬結(jié)果將更加真實和可靠，能夠更好地支持氣候預(yù)測、災(zāi)害預(yù)警和資源管理等實際應(yīng)用。第二，高性能計算技術(shù)的應(yīng)用將極大提升模擬的效率和規(guī)模。特別是隨著并行計算和分布式計算技術(shù)的不斷發(fā)展，模擬計算將能夠更好地利用多核處理器和集群系統(tǒng)的計算資源，實現(xiàn)更高的計算速度和更大的計算規(guī)模。這使得研究人員能夠更快地獲得模擬結(jié)果，并且能夠處理更大的數(shù)據(jù)規(guī)模和更復(fù)雜的模擬問題。第三，算法精度和穩(wěn)定性的提升將是未來發(fā)展的關(guān)鍵。為了適應(yīng)更加復(fù)雜的模擬需求，算法必須具備更高的精度和穩(wěn)定性。研究人員需要不斷優(yōu)化和改進(jìn)算法，提高其計算效率和結(jié)果的準(zhǔn)確性。第四，網(wǎng)格優(yōu)化和自適應(yīng)算法的研究將是提高計算效率和結(jié)果準(zhǔn)確性的重要手段。通過優(yōu)化網(wǎng)格結(jié)構(gòu)和自適應(yīng)算法，模擬計算能夠更好地適應(yīng)復(fù)雜的地形和邊界條件，提高計算結(jié)果的準(zhǔn)確性和可靠性。第五，數(shù)據(jù)驅(qū)動和機(jī)器學(xué)習(xí)技術(shù)的融合將為地球動力學(xué)研究帶來新的視角和方法。通過將機(jī)器學(xué)習(xí)算法應(yīng)用于模擬數(shù)據(jù)，研究人員能夠發(fā)現(xiàn)新的模式和關(guān)系，其應(yīng)用包括提高氣候預(yù)測和災(zāi)害預(yù)警的準(zhǔn)確性等密切關(guān)乎人類生活的場景。

地球動力學(xué)數(shù)值計算方法的應(yīng)用與發(fā)展是一個不斷進(jìn)化的領(lǐng)域。隨著計算技術(shù)的持續(xù)進(jìn)步和算法的不斷創(chuàng)新，我們相信數(shù)值模擬將在解決地球動力學(xué)領(lǐng)域的重大科學(xué)問題和推動該領(lǐng)域的發(fā)展中發(fā)揮更加關(guān)鍵的作用。未來的研究將更加注重跨學(xué)科的合作，以及新算法和新技術(shù)的開發(fā)，以期達(dá)到對地球甚至行星演化更深層次的理解和更精確的預(yù)測。

參考文獻(xiàn)（References）：

[1] Xiong F， Zhu C， Feng G， et al. A Three-Dimensional Coupled Thermo-Hydro Model for Geothermal Development in Discrete Fracture Networks of Hot Dry Rock Reservoirs[J]. Gondwana Research， 2023， 122： 331-347.

[2] Liu J， Yang H， Xu K， et al. Genetic Mechanism of Transfer Zones in Rift Basins： Insights from Geomechanical Models[J]. GSA Bulletin， 2022， 134（9/10）： 2436-2452.

[3] Liu Z， Zheng L， Zuo Y， et al. Investigation of Three-Dimensional Model Reconstruction and Fractal Characteristics of Crack Propagation in Jointed Sandstone[J]. Geomechanics and Geophysics for Geo-Energy and Geo-Resources， 2024， 10（1）： 75.

[4] Kulik D A， Wagner T， Dmytrieva S V， et al. GEM-Selektor Geochemical Modeling Package： Revised Algorithm and GEMS3K Numerical Kernel for Coupled Simulation Codes[J]. Computational Geosciences， 2013， 17（1）： 1-24.

[5] Steefel C I， Appelo C A J， Arora B， et al. Reactive Transport Codes for Subsurface Environmental Simulation[J]. Computational Geosciences， 2015， 19（3）： 445-478.

[6] Zhang X， Wu Y， Zhai E， et al. Coupling Analysis of the Heat-Water Dynamics and Frozen Depth in a Seasonally Frozen Zone[J]. Journal of Hydrology， 2021， 593（1）： 125603.

[7] Tsuboi S， Ando K， Miyoshi T， et al. A 1.8 Trillion Degrees-of-Freedom， 1.24 Petaflops Global Seismic Wave Simulation on the K Computer[J]. The International Journal of High Performance Computing Applications， 2016， 30（4）： 411-422.

[8] Sabah M， Ameri M J， Hofmann H， et al. Numerical Modeling of Injection-Induced Earthquakes Based on Fully Coupled Thermo-Poroelastic Boundary Element Method[J]. Geothermics， 2022， 105： 102481.

[9] Yuen R， De Vahl Davis G， Leona E. The Influence of Moisture on the Combustion of Wood[J]. Numerical Heat Transfer， Part A： Applications， 2000， 38（3）： 257-280.

[10] Gmeiner B， Rüde U， Stengel H. Performance and Scalability of Hierarchical Hybrid Multigrid Solvers for Stokes"" Systems[J]. SIAM Journal on Scientific Computing， 2015， 37（2）： C143-C168.

[11] Omlin S， Rss L， Podladchikov Y Y. Simulation of Three-Dimensional Viscoelastic Deformation Coupled to Porous Fluid Flow[J]. Tectonophysics， 2018， 746： 695-701.

[12] Rss L， Duretz T， Podladchikov Y Y. Resolving Hydromechanical Coupling in Two and Three Dimensions： Spontaneous Channelling of Porous Fluids Owing to Decompaction Weakening[J]. Geophysical Journal International， 2019， 218（3）： 1591-1616.

[13] Morra G， Yuen D A， Tufo H M， et al. Fresh Outlook on Numerical Methods for Geodynamics： Part 2： Big Data， HPC， Education[M]. Amsterdam： Elsevier， 2021.

[14] Rene G， Juliane D， Wolfgang B， et al. On Formulations of Compressible Mantle Convection[J]. Geophysical Journal International，2020， 221（2）： 1264-1280.

[15] Gerya T. Introduction to Numerical Geodynamic Modelling[M]. Zürich： Cambridge University Press， 2019.

[16] 劉春梅，鐘柳強(qiáng)，舒適，等. 平面彈性問題自適應(yīng)有限元方法的收斂性分析[J]. 應(yīng)用數(shù)學(xué)和力學(xué)，2014，35（9）：969-978.

Liu" Chunmei， Zhong Liuqiang， Shu Shi， et al. Convergence of an Adaptive Finite Element Method for 2D Elasticity Problems[J]. Applied Mathematics and Mechanics， 2014， 35（9）： 969-978.

[17] Barlow J. Optimal Stress Locations in Finite Elementmodels[J]. International Journal for Numerical Methods in Engineering， 1976（10）： 243-251.

[18] Strang G，" Fix G J. An Analysis of the Finite Element Method[J]. Applied Mathematics and Mechanics， 2010， 55（11）：696-697.

[19] Newton I S， Huygens C， Motte A， et al. Mathematical Principles of Natural Philosophy[M]. London： H D Symonds， 1803.

[20] Gander M J， Wanner G. From Euler， Ritz， and Galerkin to Modern Computing[J]. SIAM Review， 2012， 54（4）： 627-666.

[21] Hrennikoff A. Solution of Problems of Elasticity by the Framework Method[J]. Journal of Applied Mechanics， 1941， 8（4）： A169-A175.

[22] Courant R. Variational Methods for the Solution of Problems of Equilibrium and Vibrations[J]. Bulletin of the American Mathematical Society， 1943， 49（1）： 1-24.

[23] Turner M J， Clough R W， Martin H C， et al. Stiffness and Deflection Analysis of Complex Structures[J]. Journal of the Aeronautical Sciences， 1956， 23（9）： 805-823.

[24] Clough R W. The Finite Element Method in Plane Stress Analysis[C]//2nd Conference on Electronic Computation. Pittsburgh， PA： American Society of Civil Engineers， 1960： 345-378.

[25] 馮康. 基于變分原理的差分格式[J]. 應(yīng)用數(shù)學(xué)和計算數(shù)學(xué)，1965（1/2/3/4）： 238-262.

Feng Kang. Difference Schemes Based on Variational Principles[J]. Acta Mathematicae Applicatae Sinica， 1965（1/2/3/4）： 238-262.

[26] 臧紹先，寧杰遠(yuǎn). 俯沖帶的負(fù)浮力及其影響因素[J]. 地球物理學(xué)報，1994， 37（2）： 174-183.

Zang Shaoxian， Ning Jieyuan. Subduction Zone Negative Buoyancy and Its Influencing Factors[J]. Chinese Journal of Geophysics， 1994， 37（2）： 174-183.

[27] 石耀霖，張健. 活動海嶺俯沖與島弧火山活動的熱模擬研究[J]. 地球物理學(xué)報，1998， 41（2）： 174-181.

Shi Yaolin， Zhang Jian. Thermal Simulation Study on the Subduction of Active Ridges and Volcanic Activity of Island Arcs[J]. Chinese Journal of Geophysics， 1998， 41（2）： 174-181.

[28] 劉亞靜，葉國揚(yáng)，毛興華，等. 俯沖帶深部應(yīng)力場的二維粘彈性有限元數(shù)值模擬[J]. 地震學(xué)報，2002， 24（3）： 285-292.

Liu Yaijing， Ye Guoyang， Mao Xinghua， et al. 2-D Viscoelastic FEM Simulation on Stress State in the Deep Part of a Subducted Slab[J]. Earthquake Science， 2002， 24（3）： 285-292.

[29] 唐湘蓉，李晶. 構(gòu)造應(yīng)力場有限元數(shù)值模擬在裂縫預(yù)測中的應(yīng)用[J]. 特種油氣藏，2005， 12（2）： 25-28.

Tang Xiangrong， Li Jing. Application of Finite Element Numerical Simulation of Tectonic Stress Field in Fracture Prediction[J]. Special Oil amp; Gas Reservoirs， 2005， 12（2）： 25-28.

[30] 沈海超，程遠(yuǎn)方，王京印，等. 斷層對地應(yīng)力場影響的有限元研究[J]. 大慶石油地質(zhì)與開發(fā)，2007， 26（2）： 34-37.

Shen Haichao， Cheng Yuanfang， Wang Jingyin， et al. Study of Finite Element on Effects of Faults on Ground Stress Field[J]. Petroleum Geology and Oilfield Development in Daqing， 2007， 26（2）： 34-37.

[31] 白玉柱，周慶. 逆沖斷層錯動形成變形場及應(yīng)力場的有限元模型：以映秀-青川斷層為例[J]. 地震地質(zhì)，2013， 35（4）： 721-730.

Bai Yuzhu， Zhou Qing. The Study on the Finite Element Model of Deformation and Stress Fields Due to the Structural Motion of Inverse Fault： A Case of Yingxiu-Qingchuan Fault[J]. Seismology and Geology， 2013， 35（4）： 721-730.

[32] Davies D R， Davies J H， Hassan O， et al. Adaptive Finite Element Methods in Geodynamics： Convection Dominated Mid-Ocean Ridge and Subduction Zone Simulations[J]. International Journal of Numerical Methods for Heat amp; Fluid Flow， 2008， 18（7/8）： 1015.

[33] 戴黎明，李三忠，樓達(dá)，等. 亞洲大陸主要活動塊體的現(xiàn)今構(gòu)造應(yīng)力數(shù)值模擬[J]. 吉林大學(xué)學(xué)報（地球科學(xué)版）， 2013， 43（2）： 469-483.

Dai Liming， Li Sanzhong， Lou Da， et al. Numerical Modeling of Present-Day Structural Stress of Major Active Blocks in the Asian Continent[J]. Journal of Jilin University （Earth Science Edition）， 2013， 43（2）： 469-483.

[34] 袁杰，朱守彪. 斷層自發(fā)破裂動力過程的有限單元法模擬[J]. 地球物理學(xué)報，2014， 57（1）： 138-156.

Yuan Jie， Zhu Shoubiao. FEM Simulation of the Dynamic Processes of Fault Spontaneous Rupture[J]. Chinese Journal of Geophysics， 2014， 57（1）： 138-156.

[35] 王仁，何國琦，殷有泉，等. 華北地區(qū)地震遷移規(guī)律的數(shù)學(xué)模擬[J]. 地震學(xué)報，1980， 2（1）： 32-42.

Wang Ren， He Guoqi， Yin Youquan， et al. Mathematical Simulation of the Seismic Migration Pattern in the North China Region[J]. Earthquake Science， 1980， 2（1）： 32-42.

[36] 王仁，黃杰藩，孫荀英，等. 華北地震構(gòu)造應(yīng)力場的模擬[J]. 中國科學(xué)：地球科學(xué)：B輯， 1982，12 （4）： 337-344.

Wang Ren， Huang Jiefan， Sun Xunying， et al. Simulation of the Stress Field in the Seismotectonic Structure of North China[J]. Science China：" Earth Sciences： Series B， 1982， 12（4）： 337-344.

[37] 張郢珍，粟生平，梁北援. 魯南地區(qū)歷史地震活動特征及未來地震趨勢的有限元分析[J]. 中國地震， 1988， 4（3）： 92-95.

Zhang Yingzhen， Su Shengping， Liang Beiyuan. Finite Element Analysis of the Characteristics of Historical Earthquake Activity and Future Seismic Trends in the Southern Shandong Region[J]. Earthquake Research in China， 1988， 4（3）： 92-95.

[38] 梅世蓉，梁北援. 唐山地震孕震過程的數(shù)值模擬[J]. 中國地震， 1989， 5（3）： 11-19.

Mei Shirong， Liang Beiyuan. Numerical Simulation of the Seismic Rupture Process of the Tangshan Earthquake[J]. Earthquake Research in China， 1989， 5（3）： 11-19.

[39] 王繼存，黃清陽，續(xù)春榮，等. 川滇菱形塊構(gòu)造應(yīng)力場的數(shù)值模擬[J]. 地震地質(zhì)，1991， 13（1）： 67-72.

Wang Jicun， Huang Qingyang， Xu Chunrong， et al. Numerical Simulation of the Stress Field in the Sichuan-Yunnan Rhomboid Block Tectonics[J]. Seismology and Geology， 1991， 13（1）： 67-72.

[40] 梁海華，侯建軍，劉樹文，等. 中國構(gòu)造應(yīng)力場與大震復(fù)發(fā)周期關(guān)系的數(shù)值模擬[J]. 地震地質(zhì)，1999， 21（1）： 51-57.

Liang Haihua， Hou Jianjun， Liu Shuwen， et al. Numerical Simulation of the Relationship Between the Tectonic Stress Field and the Recurrence Interval of Major Earthquakes in China[J]. Seismology and Geology， 1999， 21（1）： 51-57.

[41] Parsons T. Post-1906 Stress Recovery of the San Andreas Fault System Calculated from Three-Dimensional Finite Element Analysis[J]. Journal of Geophysical Research： Solid Earth， 2002， 107（B8）： 2162.

[42] 朱守彪，邢會林，謝富仁，等. 地震發(fā)生過程的有限單元法模擬：以蘇門答臘俯沖帶上的大地震為例[J]. 地球物理學(xué)報，2008， 51（2）： 460-468.

Zhu Shoubiao， Xing Huilin， Xie Furen， et al. Simulation of Earthquake Processes by Finite Element Method： The Case of Megathrust Earthquakes on the Sumatra SubductionZone[J]. Chinese Journal of Geophysics， 2008， 51（2）： 460-468.

[43] 李平恩，廖力，奉建州，等. 1900年以來巴顏喀拉塊體應(yīng)力演化與周緣強(qiáng)震關(guān)系的數(shù)值模擬研究[J]. 地球物理學(xué)報，2019， 62（11）： 4170-4188.

Li Ping’en， Liao Li， Feng Jianzhou， et al. Numerical Simulation of Relationship Between Stress Evolution and Strong Earthquakes Around the Bayan Har Block Since 1900[J]. Chinese Journal of Geophysics， 2019， 62（11）： 4170-4188.

[44] Qu W， Gao Y， Zhang Q， et al. Present Crustal Deformation and Stress-Strain Fields of North China Revealed from GPS Observations and Finite Element Modelling[J]. Journal of Asian Earth Sciences， 2019， 183（Oct.1）： 103959.

[45] Bahuguna A， Shanker D. Simulation of Intraplate Stress Distribution of the Indian Tectonic Plate Using the Finite Element Method[J].Pure and Applied Geophysics， 2022， 179（1）：125-148.

[46] Karabulut M F， Gülal V E. Accuracy Assessment of Mesh Types in Tectonic Modeling to Reveal Seismic Potential Using the Finite Element Method： A Case Study at the Main Marmara Fault[J]. Applied Sciences， 2023， 13（24）： 13297.

[47] 馬峰，高俊，王貴玲，等. 雄安新區(qū)容城地?zé)崽锾妓猁}巖熱儲采灌數(shù)值模擬[J]. 吉林大學(xué)學(xué)報（地球科學(xué)版）， 2023， 53（5）： 1534-1548.

Ma Feng， Gao Jun， Wang Guiling， et al. Numerical Simulation of Exploitation and Reinjection of Carbonate Geothermal Reservoir in Rongcheng Geothermal Field， Xiong’an New Area[J]. Journal of Jilin University （Earth Science Edition）， 2023， 53（5）： 1534-1548.

[48] Xu F， Hajibeygi H， Sluys L J. Adaptive Multiscale Extended Finite Element Method （MS-XFEM） for the Simulation of Multiple Fractures Propagation in Geological Formations[J]. Journal of Computational Physics， 2023， 486（C）： 112-114.

[49] Chen Y， Yang H， Ye Y， et al. Generation of 3D Finite Element Mesh of Layered Geological Bodies in Intersecting Fault Zones[J]. Plos One， 2024， 19（1）： e0293193.

[50] Euler" L. Institutionum Calculi Integralis Volumen Primum[M]. St Petersburg： Imperial Academy of Sciences， 1768.

[51] Runge C， Precht. ber die Beziehungen des Atomgewichts der Elemente zu Deren Spektrum[J]. Zeitschrift für Analytische Chemie， 1905， 44（3）： 265.

[52] Alterman Z S， Karal F C J. Propagation of Elastic Waves in Layered Media by Finite Difference Methods[J]. Bulletin of the Seismological Society of America， 1968， 58（1）： 367-398.

[53] Virieux J. P-SV Wave Propagation in Heterogeneous Media： Velocity-Stress Finite-Difference Method[J]. Geophysics， 1986， 51（4）： 889-901.

[54] Dablain M A. The Application of High-Order Differencing to the Scalar Wave Equation[J]. Geophysics， 1986， 51（1）： 54-66.

[55] Bayliss A， Jordan K E， LeMesurier B J， et al. A Fourth-Order Accurate Finite-Difference Scheme for the Computation of Elastic Waves[J]. Bulletin of the Seismological Society of America， 1986， 76（4）： 1115-1132.

[56] Levander A R. Fourth-Order Finite-Difference P-SV Seismograms[J]. Geophysics， 1988， 53（11）： 1425-1436.

[57] Jastram C， Tessmer E. Elastic Modelling on a Grid with Vertically Varying Spacing[J]. Geophysical Prospecting， 1994， 42（4）： 357-370.

[58] Falk J， Tessmer E， Gajewski D. Efficient Finite-Difference Modelling of Seismic Waves Using Locally Adjustable Time Steps[J]. Geophysical Prospecting， 1998， 46（6）： 603-616.

[59] Graves R W. Simulating Seismic Wave Propagation in 3D Elastic Media Using Staggered-Grid Finite Differences[J]. Bulletin of the Seismological Society of America， 1996， 86（4）： 1091-1106.

[60] Boore D. Finite Difference Methods for Seismic Wave Propagation in Heterogeneous Materials[M]. Amsterdam： Elsevier， 1972.

[61] Kelly K R， Ward R W， Treitel S， et al. Synthetic Seismograms： A Finite-Difference Approach[J]. Geophysics， 1976， 41（1）： 2-27.

[62] Madariaga R. Dynamics of an Expanding Circular Fault[J]. Bulletin of the Seismological Society of America， 1976， 66（3）： 639-666.

[63] 高孟潭，俞言祥，張曉梅，等. 北京地區(qū)地震動的三維有限差分模擬[J]. 中國地震，2002， 18（4）： 356-364.

Gao Mengtan， Yu Yanxiang， Zhang Xiaomei， et al. 3D Finite-Difference Simulation of Ground Motion in the Beijing Area[J]. Earthquake Research in China， 2002， 18（4）： 356-364.

[64] Hernlund J W， Tackley P J. Three-Dimensional Spherical Shell Convection at Infinite Prandtl Number Using the Cubed Sphere' Method[M]. Amsterdam： Elsevier， 2003.

[65] Igel H， Riollet B， Mora P. Accuracy of Staggered 3-D Finite-Difference Grids for Anisotropic Wave Propagation[C]//SEG Technical Program Expanded Abstracts 1992. Houston： Society of Exploration Geophysicists， 1992： 1244-1246.

[66] Pitarka A. 3D Elastic Finite-Difference Modeling of Seismic Motion Using Staggered Grids with Nonuniform Spacing[J]. Bulletin of the Seismological Society of America， 1999， 89（1）： 54-68.

[67] Bohlen T， Saenger E H. Accuracy of Heterogeneous Staggered-Grid Finite-Difference Modeling of Rayleigh Waves[J]. Geophysics， 2006， 71（4）： T109-T115.

[68] 李建平，翁愛華，李世文，等. 基于球坐標(biāo)系下有限差分的地磁測深三維正演[J]. 吉林大學(xué)學(xué)報（地球科學(xué)版）， 2018， 48（2）： 411-419.

Li Jianping， Weng Aihua， Li Shiwen， et al. 3-D Forward Method for Geomagnetic Depth Sounding Based on Finite Difference Method in Spherical Coordinate[J]. Journal of Jilin University （Earth Science Edition）， 2018， 48（2）： 411-419.

[69] 韓麗，胥鐵瀟，金勝昔，等. 基于雙電層模型的頻域震電響應(yīng)數(shù)值模擬[J]. 吉林大學(xué)學(xué)報（地球科學(xué)版）， 2022， 52（3）： 737-743.

Han Li， Xu Tiexiao， Jin Shengxi， et al. Numerical Simulation of Seismoelectric Response inthe Frequency Domain Based on Double Electric Layer Model[J]. Journal of Jilin University （Earth Science Edition）， 2022， 52（3）： 737-743.

[70] 張志佳，孫章慶，王瑞湖，等. 復(fù)雜地表起伏多重變加密網(wǎng)格有限差分法波動方程數(shù)值模擬[J]. 吉林大學(xué)學(xué)報（地球科學(xué)版）， 2023， 53（2）： 598-608.

Zhang Zhijia， Sun Zhangqing， Wang Ruihu， et al. Numerical Simulation of Wave Equation Using Finite Difference Method with Rugged Variable Refined Multigrid in Complex Surface[J]. Journal of Jilin University （Earth Science Edition）， 2023， 53（2）： 598-608.

[71] Gerya T V， May D A， Duretz T. An Adaptive Staggered Grid Finite Difference Method for Modeling Geodynamic Stokes Flows with Strongly Variable Viscosity[J]. Geochemistry， Geophysics， Geosystems， 2013， 14（4）： 1200-1225.

[72] Chávez-Negrete C， Domínguez-Mota F J， Román-Gutiérrez R. Interface Formulation for Generalized Finite Difference Method for Solving Groundwater Flow[J]. Computers and Geotechnics， 2024， 166： 105990.

[73] Jaswon M A. Integral Equation Methods in Potential Theory： I[J]. Proceedings of the Royal Society A： Mathematical， Physical and Engineering Sciences， 1963， 275： 23-32.

[74] Rizzo F J. An Integral Equation Approach to Boundary Value Problems of Classical Elastostatics[J]. Quarterly of Applied Mathematics， 1967， 25（1）： 83-95.

[75] Sweedlow J L， Cruse T A. Formulation of Boundary Integral Equations for Three-Dimensional Elasto-Plastic Flow[J]. International Journal of Solids and Structures， 1971， 7（12）： 1673-1683.

[76] Cruse T A， Vanburen W. Three-Dimensional Elastic Stress Analysis of a Fracture Specimen with an Edge Crack[J]. International Journal of Fracture Mechanics， 1971， 7（1）： 1-15.

[77] Mendelson A. Solution of Elastoplastic Torsion Problem by Boundary Integral Method [R]. Washington D C： NASA， 1975.

[78] Yin T， Hsiao G C， Xu L. Boundary Integral Equation Methods for the Two-Dimensional Fluid-Solid Interaction Problem[J]. SIAM Journal on Numerical Analysis， 2017， 55（5）： 2361-2393.

[79] Quevedo L， Morra G， Müller R D. Parallel Fast Multipole Boundary Element Method for Crustal Dynamics[J]. IOP Conference Series： Materials Science and Engineering， 2010， 10（1）： 012012.

[80] Liu Y， Wen Y， Li Z， et al. Coseismic Fault Model of the 2017 M_W 6.5 Jiuzhagou Earthquake and Implications for the Regional Fault Slip Pattern[J]. Geodesy and Geodynamics， 2022， 13（2）： 104-113.

[81] 杜慶華，姚振漢，岑章志.我國工程中邊界元法研究的十年[J].力學(xué)與實踐， 1989， 11（1）：5.

Du Qinghua， Yao Zhenhan， Cen Zhangzhi. A Decade of Boundary Element Method Research in Chinese Engineering [J]. Mechanics in Engineering， 1989， 11（1）： 5.

[82] 董春迎. 彈塑性邊界元法的若干基礎(chǔ)性研究及在接觸問題上的應(yīng)用[D]. 北京：清華大學(xué)，1992.

Dong Chunying. Some Basic Studies on Elastoplastic Boundary Element Method and Its Application in Contact Problems [D]. Beijing： Tsinghua University， 1992.

[83] 金朝嵩. 自適應(yīng)邊界元法的后驗誤差估計[J]. 重慶建筑大學(xué)學(xué)報，2000， 22（6）： 16-19.

Jin Chaosong. A Posteriori Error Estimates for Adaptive Boundary Element Methods[J]. Journal of Chongqing University， 2000， 22（6）： 16-19.

[84] 吳云，申重陽，周碩愚，等. 基于邊界元的非連續(xù)（塊體系統(tǒng)）形變反分析法[J]. 武漢大學(xué)學(xué)報（信息科學(xué)版）， 2003， 28（3）： 345-350.

Wu Yun， Shen Chongyang， Zhou Shuoyu， et al. An Inversion Method of DDA with BEM[J]. Geomatics and Information Science of Wuhan University， 2003， 28（3）： 345-350.

[85] 焦龍，吳云，張躍剛. 用邊界元法研究川西地區(qū)的形變場與應(yīng)變場[J]. 大地測量與地球動力學(xué)，2010， 30（3）： 35-38.

Jiao Long， Wu Yun， Zhang Yuegang. Research on Deformation Field and Strain Field in Western Sichuan Area by BEM[J]. Geodesy and Geodynamics， 2010， 30（3）： 35-38.

[86] Ogawa M， Schubert G， Zebib A. Numerical Simulations of Three-Dimensional Thermal Convection in a Fluid with Strongly Temperature-Dependent Viscosity[J]. Journal of Fluid Mechanics， 1991， 233（1）： 299-328.

[87] Tackley P J. Self-Consistent Generation of Tectonic Plates in Time-Dependent， Three-Dimensional Mantle Convection Simulations[J]. Geochemistry， Geophysics， Geosystems， 2000， 1（8）： 1021.

[88] Gerya T V， Yuen D A. Characteristics-Based Marker-in-Cell Method with Conservative Finite-Differences Schemes for Modeling Geological Flows with Strongly Variable Transport Properties[J]. Physics of the Earth and Planetary Interiors， 2003， 140（4）： 293-318.

[89] Choblet G， [KG-0.5mm]C[DD（-1][HT5”]ˇ[DD）]adek O， Couturier F， et al. DIPUS： A New Tool to Study the Dynamics of Planetary Interiors[J]. Geophysical Journal International， 2007， 170（1）： 9-30.

[90] Aziz K， Settari A. Petroleum Reservoir Simulation[M]. Radarweg： Elsevier Applied Science Publishers， 1986.

[91] Jung B， Garven G， Boles J R. The Geodynamics of Faults and Petroleum Migration in the Los Angeles Basin， California[J]. American Journal of Science， 2015， 315（5）： 412-459.

[92] 馬瑞杰，李欣. 巖溶裂隙地下水流數(shù)學(xué)模型求解的有限體積法及應(yīng)用[J]. 吉林大學(xué)學(xué)報（地球科學(xué)版）， 2005， 35（6）： 84-87.

Ma Ruijie， Li Xin. The Finite Volume Method and Application of Mathematical Model of Karst Groundwater Flow[J]. Journal of Jilin University （Earth Science Edition）， 2005， 35（6）： 84-87.

[93] Jenny P， Lee S H， Tchelepi H A. Adaptive Fully Implicit Multi-Scale Finite-Volume Method for Multi-Phase Flow and Transport in Heterogeneous Porous Media[J]. Journal of Computational Physics， 2006， 217（2）： 627-641.

[94] Zhang Z. Error Estimates of Finite Volume Element Method for the Pollution in Groundwater Flow[J]. Numerical Methods for Partial Differential Equations， 2009， 25（2）： 259-274.

[95] Cui X， Wong L N Y. A 3D Thermo-Hydro-Mechanical Coupling Model for Enhanced Geothermal Systems[J]. International Journal of Rock Mechanics and Mining Sciences， 2021， 143（6）： 104744.

[96] Huang Y， Zhou Z F， Yu Z B. Simulation of Solute Transport Using a Coupling Model Based on Finite Volume Method in Fractured Rocks[J]. Journal of Hydrodynamics， 2010， 22（1）： 129-136.

[97] Lucy L B. A Numerical Approach to the Testing of the Fission Hypothesis[J]. The Astronomical Journal， 1977， 82（12）： 1013-1024.

[98] Nayroles B， Touzot G， Villon P. Generalizing the Finite Element Method： Diffuse Approximation and Diffuse Elements[J]. Computational Mechanics， 1992， 10（5）： 307-318.

[99] Belytschko T， Lu Y Y， Gu L. Element-Free Galerkin Methods[J]. International Journal for Numerical Methods in Engineering， 1994， 37（2）： 229-256.

[100] Liu W K， Jun S， Zhang Y F. Reproducing Kernel Particle Methods[J]. International Journal for Numerical Methods in Fluids， 1995， 20（8）： 1081-1106.

[101] Hardy R L. Multiquadric Equations of Topography and Other Irregular Surfaces[J]. Journal of Geophysical Research， 1971， 76（8）： 1905-1915.

[102] 張雄，劉巖，馬上. 無網(wǎng)格法的理論及應(yīng)用[J]. 力學(xué)進(jìn)展，2009， 39（1）： 1-36.

Zhang Xiong， Liu Yan， Ma Shang. The Theory and Application of Meshfree Methods[J]. Progress in Mechanics， 2009， 39（1）： 1-36.

[103] Nguyen V P， Rabczuk T， Bordas S， et al. Meshless Methods： A Review and Computer Implementation Aspects[J]. Mathematics and Computers in Simulation， 2008， 79（3）： 763-813.

[104] 焦玉玲，劉金英，楊天行，等. 無網(wǎng)格法在地下水水位預(yù)測中的應(yīng)用[J]. 吉林大學(xué)學(xué)報（地球科學(xué)版）， 2007， 37（3）： 535-540.

Jiao Yuling， Liu Jinying， Yang Tianxing， et al. Application of Meshless Method in Groundwater Level Forecast[J]. Journal of Jilin University （Earth Science Edition）， 2007， 37（3）： 535-540.

[105] 王月英. 地震波正演模擬中無單元Galerkin法初探[J]. 地球物理學(xué)進(jìn)展，2007， 22（5）： 1539-1544.

Wang Yueying. Study of Element-Free Galerkin Method in Seismic Forward Modeling[J]. Progress in Geophysics， 2007， 22（5）： 1539-1544.

[106] 李坤，黃其柏，林立廣. 二階時域波動方程的無網(wǎng)格方法求解[J]. 華中科技大學(xué)學(xué)報（自然科學(xué)版）， 2011， 39（3）： 26-29.

Li Kun， Huang Qibai， Lin Liguang. Solving Second-Order Time Domain Wave Equations by Meshless Method[J]. Journal of Huazhong University of Science and Technology （Natural Science Edition）， 2011， 39（3）： 26-29.

[107] Wang L， Wang Z， Qian Z. A Meshfree Method for Inverse Wave Propagation Using Collocation and Radial Basis Functions[J]. Computer Methods in Applied Mechanics and Engineering， 2017， 322： 311-350.

[108] 馮德山，王洪華，戴前偉. 基于無單元Galerkin法探地雷達(dá)正演模擬[J]. 地球物理學(xué)報，2013， 56（1）： 298-308.

Feng Deshan， Wang Honghua， Dai Qianwei. Forward Simulation of Ground Penetrating Radar Based on the Element-Free Galerkin Method[J]. Chinese Journal of Geophysics， 2013， 56 （1）： 298-308.

[109] 嵇艷鞠，黃廷哲，黃婉玉，等. 起伏地形下各向異性的2D大地電磁無網(wǎng)格法數(shù)值模擬[J]. 地球物理學(xué)報，2016， 59（12）： 4483-4493.

Ji Yanju， Huang Tingzhe， Huang Wanyu， et al. 2D Anisotropic Magnetotelluric Numerical Simulation Using Meshfree Method Under Undulating Terrain[J]. Chinese Journal of Geophysics， 2016， 59（12）： 4483-4493.

[110] 李俊杰，嚴(yán)家斌. 無網(wǎng)格點插值法大地電磁二維正演數(shù)值模擬[J]. 石油物探，2014， 53（5）： 617-626.

Li Junjie， Yan Jiabin. Magnetotelluric Two-Dimensional Forward Numerical Modeling by Meshfree Point Interpolation Method[J]. Geophysical Prospecting for Petroleum， 2014， 53（5）： 617-626.

[111] 魏德敏，胡源喆，江雪玲. 邊坡彈塑性穩(wěn)定性分析的無網(wǎng)格法[J]. 力學(xué)與實踐，2012， 34（3）： 12-17.

Wei Demin， Hu Yuanzhe， Jiang Xueling. Application of Meshless Method in Analysis of Elastic-Plastic Stability of an Edge Slope[J]. Mechanics in Engineering， 2012， 34（3）： 12-17.

[112] McDougall S， Hungr O. A Model for the Analysis of Rapid Landslide Motion Across Three-Dimensional Terrain[J]. Canadian Geotechnical Journal， 2004， 41（6）： 1084-1097.

[113] Lin S， Cao X， Zheng H， et al. An Improved Meshless Numerical Manifold Method for Simulating Complex Boundary Seepage Problems[J]. Computers and Geotechnics， 2023， 155： 105211.

[114] Cao X， Lin S， Liang Z， et al. Meshless Numerical Manifold Method with Novel Subspace Tracking and CSS Locating Techniques for Slope Stability Analysis[J]. Computers and Geotechnics， 2024， 166： 106025.

[115] Van Keken P E， Wada I， Sime N， et al. Thermal Structure of the Forearc in Subduction Zones： A Comparison of Methodologies[J]. Geochemistry， Geophysics， Geosystems， 2019， 20（7）： 3268-3288.

[116] Aagaard B T， Knepley M G， Williams C A. A Domain Decomposition Approach to Implementing Fault Slip in Finite-Element Models of Quasi-Static and Dynamic Crustal Deformation[J]. Journal of Geophysical Research： Solid Earth， 2013， 118（6）： 3059-3079.

[117] Ward L A， Smith-Konter B R， Xu X， et al. Seismic Moment Accumulation Response to Lateral Crustal Variations of the San Andreas Fault System[J]. Journal of Geophysical Research： Solid Earth， 2021， 126（4）： e2020JB021208.

[118] Holden L， Cas R， Fournier N， et al. Modelling Ground Deformation Patterns Associated with Volcanic Processes at the Okataina Volcanic Centre[J]. Journal of Volcanology and Geothermal Research， 2017， 344： 65-78.

[119] Aochi H， Tsuda K. Dynamic Rupture Simulations Based on Depth-Dependent Stress Accumulation[J]. Geophysical Journal International， 2022， 233（1）： 182-194.

[120] Oral E， Ampuero J P， Ruiz J， et al. A Method to Generate Initial Fault Stresses for Physics-Based Ground Motion Prediction Consistent with Regional Seismicity[J]. Bulletin of the Seismological Society of America， 2022， 112（6）： 2812–2827.

[121] 祝賀君，劉沁雅，楊繼東. 地震學(xué)全波形反演進(jìn)展[J]. 地球與行星物理論評， 2023， 54（3）： 287-317.

Zhu Hejun， Liu Qinya， Yang Jidong. Recent Progress on Full Waveform Inversion[J]. Reviews of Geophysics and Planetary Physics， 2023， 54（3）： 287-317.

[122] 陳衛(wèi)忠，伍國軍，戴永浩，等. 廢棄鹽穴地下儲氣庫穩(wěn)定性研究[J]. 巖石力學(xué)與工程學(xué)報，2006， 25（4）： 848-854.

Chen Weizhong， Wu Guojun， Dai Yonghao， et al. Stability Analysis of Abandoned Salt Caverns Used for Underground Gas Storage[J]. Chinese Journal of Rock Mechanics and Engineering， 2006， 25（4）： 848-854.

[123] 張廣明，熊春明，劉合，等. 復(fù)雜斷塊地應(yīng)力場數(shù)值模擬方法研究[J]. 斷塊油氣田，2011， 18（6）： 710-713.

Zhang Guangming， Xiong Chunming， Liu He， et al. Numerical Simulation Method for In-Situ Stress Field in Complex Fault Block[J]. Fault-Block Oil amp; Gas Field， 2011， 18（6）： 710-713.

[124] 李帆，楊建國. 黃土邊坡穩(wěn)定性分析方法研究[J]. 鐵道工程學(xué)報，2008， 30（12）： 33-36.

Li Fan， Yang Jianguo. Study of Analysis Method for Loess Slope"" Stability[J]. Journal of Railway Engineering Society， 2008， 30（12）： 33-36.

[125] Hampel A， Lüke J， Krause T， et al. Finite-Element Modelling of Glacial Isostatic Adjustment （GIA）： Use of Elastic Foundations at Material Boundaries Versus the Geometrically Non-Linear Formulation[J]. Computers amp; Geosciences， 2019， 122： 1-14.

[126] Volpini C， Douglas J， Nielsen A H. Guidance on Conducting 2D Linear Viscoelastic Site Response Analysis Using a Finite Element Code[J]. Journal of Earthquake Engineering， 2021， 25（6）： 1153-1170.

[127] Zhou S， Zhuang X， Rabczuk T. A Phase-Field Modeling Approach of Fracture Propagation in Poroelastic Media[J]. Engineering Geology， 2018， 240： 189-203.

[128] Xu J， Lan W， Ren C， et al. Modeling of Coupled Transfer of Water， Heat and Solute in Saline Loess Considering Sodium Sulfate Crystallization[J]. Cold Regions Science and Technology， 2021， 189： 103335.

[129] 葉海林，鄭穎人，杜修力，等. 邊坡動力破壞特征的振動臺模型試驗與數(shù)值分析[J]. 土木工程學(xué)報，2012， 45（9）： 128-135.

Ye Hailin， Zheng Yingren， Du Xiuli， et al. Shaking Table Model Test and Numerical Analysis on Dynamic Failure Characteristics of Slope[J]. China Civil Engineering Journal， 2012， 45（9）： 128-135.

[130] Mitchell M A， Peacock J R， Burgess S D. Imaging the Magmatic Plumbing of the Clear Lake Volcanic Field Using 3-D Gravity Inversions[J]. Journal of Volcanology and Geothermal Research， 2023， 435： 107758.

[131] Heagy L J， Cockett R， Kang S， et al. A Framework for Simulation and Inversion in Electromagnetics[J]. Computers amp; Geosciences， 2017， 107： 1-19.

[132] Herrera M T， Crempien J G F， Cembrano J， et al. Seismic Cycle Controlled by Subduction Geometry： Novel 3-D Quasi-Dynamic Model of Central Chile Megathrust[J]. Geophysical Journal International， 2024， 237（2）： 772-787.

[133] Hassard P. Comparison of Lattice Boltzmann and Boundary Element Methods for Stokes and Visco-Inertial Flow in a Two-Dimensional Porous Medium[J]. Transport in Porous Media， 2023， 150（3）： 675-707.

[134] Dhanya J， Raghukanth S T G. Implication of Source Models on Tsunami Wave Simulations for 2004 （Mw 9.2） Sumatra Earthquake[J]. Natural Hazards， 2020， 104（1）： 279-304.

[135] Fisher T L， Harris R A. Reconstruction of 1852 Banda Arc Megathrust Earthquake and Tsunami[J]. Natural Hazards， 2016， 83（1）： 667-689.

[136] Guo Z， Rüpke L， Tao C. Hydrothermal Foam v1.0： A 3-D Hydrothermal Transport Model for Natural Submarine Hydrothermal Systems[J]. Geoscientific Model Development， 2020， 13（12）： 6547-6565.

[137] Douglas S， Nistor I. On the Effect of Bed Condition on the Development of Tsunami-Induced Loading on Structures Using OpenFOAM[J]. Natural Hazards， 2015， 76（1）： 1335-1356.

[138] Li Y， Qi J. Salt-Related Contractional Structure and Its Main Controlling Factors of Kelasu Structural Zone in Kuqa Depression： Insights from Physical and Numerical Experiments[J]. Procedia Engineering， 2012， 31： 863-867.

[139] Cruz L， Malinski J， Wilson A， et al. Erosional Control of the Kinematics and Geometry of Fold-and-Thrust Belts Imaged in a Physical and Numerical Sandbox[J]. Journal of Geophysical Research， 2010， 115（B9）： B09404.

[140] Sandiford D， Moresi L. Improving Subduction Interface Implementation in Dynamic Numerical Models[J]. Solid Earth， 2019， 10（3）： 969-985.

[141] Liu H， Ban S， Bédard K， et al. Characteristics of Precambrian Basement Intruded by Cretaceous Geological Intrusions in Monteregian Igneous Province and Their Impacts on Regional Thermal Structure[J]. Advances in Geo-Energy Research， 2022， 6（3）： 206-220.

[142] Li M. The Influence of Uncertain Mantle Density and Viscosity Structures on the Calculations of Deep Mantle Flow and Lateral Motion of Plumes[J]. Geophysical Journal International， 2023， 233（3）： 1916-1937.

[143] Yuan Q， Li M， Desch S J， et al. Moon-Forming Impactor as a Source of Earth’s Basal Mantle Anomalies[J]. Nature， 2023， 623： 95-99.

[144] Xu H， Li J， Wang L， et al. Influence of Mantle Plume on Continental Rift Evolution： A Case Study of the East African Rift System[J/OL]. Petroleum Research， 2024. [2024-05-10].doi：10.1016/j.ptlrs.2024.02.001.

[145] Dong M， Hao T， Xu L， et al. Subduction Without Volcanic Arc Magma： Insights from Two Young Subduction Zones in the Western Pacific[J]. Tectonophysics， 2024， 874： 230231.

[146] Chambers E L， Bonadio R， Fullea J， et al. Determining Subsurface Temperature amp; Lithospheric Structure from Joint Geophysical-Petrological Inversion： A Case Study from Ireland[J]. Tectonophysics， 2023， 869： 230094.

[147] Mousavi N， Tatar M， Moghadam H S， et al. Integrated Geophysical-Petrological Model of Damavand Volcano， North Iran： Compositional Structure of Crust and Upper Mantle from Seismic Temperature Profiling and Gravity Anomaly Fitting[J]. Journal of Volcanology and Geothermal Research， 2023， 442： 107913.

[148] Oualid M， Ebong E D. Lithospheric and Asthenospheric Properties of the Saharan Platform Inferred from Potential Field， Geoid， and Heat Flow Data[J]. Journal of African Earth Sciences， 2024， 209（1）： 105124.

[149] Hardy S. Discrete Element Modelling of Extensional， Growth， Fault-Propagation Folds[J]. Basin Research， 2019， 31（3）： 584-599.

[150] Meng W J， Chen Z A， Bai W M. Numerical Simulation on Process of the Plume-Lithosphere Interaction[J]. Chinese Journal of Geophysics， 2015， 58（2）： 495-503.

[151] Fan B， Chen Z A， Bai W M. Numerical Simulation on the Melting Process of the Mantle Plume-Lithosphere Interaction[J]. Chinese Journal of Geophysics， 2018， 61（4）： 1341-1351.

[152] Hernández P A， Calvari S， Ramos A， et al. Magma Emission Rates from Shallow Submarine Eruptions Using Airborne Thermal Imaging[J]. Remote Sensing of Environment， 2014， 154： 219-225.

[153] Lichtenberg T， Schaefer L K， Nakajima M， et al. Geophysical Evolution During Rocky Planet Formation[J/OL]. arXiv Preprint， 2022. [2024-05-10].doi：10.48550/arXiv.2203.10023.

[154] Lichtenberg T， Krijt S. System-Level Fractionation of Carbon from Disk and Planetesimal Processing[J]. The Astrophysical Journal Letters， 2021， 913（2）： L20.

[155] Giunchi C， Sabadini R， Boschi E， et al. Dynamic Models of Subduction： Geophysical and Geological Evidence in the Tyrrhenian Sea[J]. Geophysical Journal International， 1996， 126（2）： 555-578.

[156] Song R， Cui M， Liu J， et al. A Pore-Scale Simulation on Thermal-Hydromechanical Coupling Mechanism of Rock[J]. Geofluids， 2017， 2017（1）： 1-12.

[157] Kaus B J P， Popov A， May D A. Recent Progress in Modelling 3D Lithospheric Deformation[C]//EGU General Assembly Conference Abstracts. Vienna： Copernicus Publications， 2012： 14116.

[158] Thielmann M， May D A， Kaus B J P. Discretization Errors in the Hybrid Finite Element Particle-in-Cell Method[J]. Pure and Applied Geophysics， 2014， 171（9）： 2165-2184.

[159] Adamuszek M， Dabrowski M， Schmalholz S M， et al. Estimating Rheological Parameters of Anhydrite from Folded Evaporite Sequences： Implications for Internal Dynamics of Salt Structure[C]//Geophysical Research Abstracts. Vienna： EGU General Assembly Conference， 2015： 9954.

[160] Okuwaki R， Yagi Y， Taymaz T， et al. Multi-Scale Rupture Growth with Alternating Directions in a Complex Fault Network During the 2023 South-Eastern Türkiye and Syria Earthquake Doublet[J]. Geophysical Research Letters， 2023， 50（12）： e2023GL103480.

[161] Adams A C， Stegman D R， Mohammadzadeh H， et al. Plume-Induced Delamination Initiated at Rift Zones on Venus[J]. Journal of Geophysical Research： Planets， 2023， 128（10）： e2023JE007879.

[162] Gülcher A J P， Golabek G J， Thielmann M， et al. Narrow， Fast， and “Cool” Mantle Plumes Caused by Strain-Weakening Rheology in Earth’s Lower Mantle[J]. Geochemistry， Geophysics， Geosystems， 2022， 23（10）： e2021GC010314.